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Frequency-domain filters for time-windowed gravitational waves from inspiralling compact bina- 
ries are constructed which combine the excellent performance of our previously developed time- 
domain P-approximants with the analytic convenience of the stationary phase approximation with- 
out a serious loss in event rate. These Fourier-domain representations incorporate the "edge oscil- 
lations" due to the (assumed) abrupt shut-off of the time-domain signal caused by the relativistic 
plunge at the last stable orbit. These new analytic approximations, the SPP-approximants, are 
not only effectual for detection and faithful for parameter estimation, but are also computation- 
^"v ally inexpensive to generate (and are faster by factors up to 10, as compared to the corresponding 

time-domain templates). The SPP approximants should provide data analysts the Fourier-domain 
templates for massive black hole binaries of total mass m < 40Mq , the most likely sources for LIGO 
and VIRGO. 
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^ . I. INTRODUCTION AND SUMMARY 

o 

The discovery of the first binary pulsar in 1974 |p| has had a very important impact on gravitational wave re- 
search. First, it proved the reality of gravitational radiation by measuring the orbital period decay Q entailed by 
the propagation at the velocity of light of the gravitational interaction between the two neutron stars making up the 
system Second, it provided the first experimental evidence that General Relativity correctly describes gravity 
in the strong-field regime Q. Third, it led to a shift in perception regarding the most promising sources for future 
O" 1 gravitational wave (GW) detectors, away from the then assumed, violent — but less predictable — gravitational 
?_i collapse associated with supcrnovac, to the more predictable, final inspiralling phase of compact binaries of neutron 
stars and black holes driven by gravitational radiation-reaction. This also led to the thrust in the laser interferometric 
gravitational wave detectors which are inherently broad-band rather than in the narrow-band bar detectors. 

•i-H . 

X: 

A. Data analysis algorithms for inspiral wave searches 

Consider a compact binary system like the binary pulsar after it has been inspiralling inwards for three hundred 
million years due to gravitational radiation-reaction. The inspiral waveform enters the detector bandwidth during 
the last few minutes of evolution of the binary. Our ability, in principle, to compute the waveform very accurately, 
allows us to track the gravitational wave phase and enhance the signal-to-noise ratio by integrating the signal for the 
interval during which it lasts in the detector band. This, in turn, requires a template with which the detector output 
may be filtered. Though template waveforms should, optimally, be exact copies of the expected signal, in practice, 
they are constructed by some approximation scheme and will differ from the actual signal in the detector output. 
Consequently, the overlap of template and signal waveforms will be less than if they had exactly matched, leading to 
a loss of potential events. Data analysis issues like these for inspiralling compact binaries of neutron stars and black 
holes have been formulated and addressed for the last twelve years even though interferometric gravitational 

wave detectors like the GEO600 § or LIGO § and VIRGO § are a year or three in the future. Much of the work in 
this area has addressed practical issues of direct relevance to data analysis strategies. These include: construction of 
templates for detection [|10|| , the number of templates, their placement, spacing, the required computing power and the 



storage or memory requirement 11 1, the order of post-Newtonian (PN) approximation adequate for detection [|12|-[14| 
parameter estimation by covariance matrix JT^-[l^] and Monte Carlo simulations Jl8| , determination of cosmological 



parameters |19| ], tests of general relativity ]2(|, one step versus hierarchical searches 21 1, effects of precession |2£ 
and of eccentricity |lj],|23|]. For the time-domain waveform, all of these works use the restricted post-Newtonian 
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approximation to quasi-circular inspiral. This keeps the crucial phase information to the best order of approximation 
then available 0], but restricts the amplitude to be Newtonian and the harmonic to the second harmonic of the 
orbital frequency. Such an approximation should be adequate for the on-line search of gravitational wave signals [§5] . 
Evidently, it is assumed that the offline analysis of the data will use the best available (unrestricted post-Newtonian) 
representation of the inspiral signals. 

B. Modelling inspiral waveforms 

The post- Newtonian approximation is basically a Taylor expansion (in powers of v/c) and all the above treatments 
use as building blocks the straightforward Taylor expansions in v/c of some intermediate quantities (orbital energy 
and gravitational- wave flux) . We shall refer to the templates based on such straightforward PN expansions as "Taylor 
approximants" (or simply T-approximants). The very slow convergence and oscillatory behaviour of the PN expansion, 
and therefore of the sequence of Taylor approximants, made imperative a search for better approximants for phasing. 
This prompted us |l3f (later referred to as D1S) to propose new approximants, with much improved convergence 
propertieSjfor application to gravitational-wave data analysis problems. 

In DIS 131, we showed how to construct a new type of time-domain approximant, called "P-approximants" , which 
not only converged faster and more monotonically, but were also more effectual (larger overlaps for detection) and 
faithful (smaller biases for parameter estimation) than the standard T-approximants. Our construction was two- 
pronged: on the one hand, it introduced new basic energy and flux functions, and on the other hand, it made a 
systematic use of Pade techniques (a well-known convergence-acceleration technique) to construct successive approx- 
imants of our new basic energy and flux functions. These new functions form a pivotal aspect of our construction 
and successfully handle issues related to appearance of non-rational functions in the energy function and logarithmic 
terms in the flux function that for long proved to be hurdles to the application of well-known Pade techniques to this 
problem. For initial LIGO, the 2.5PN P- approximants are likely to provide overlaps in excess of 96.5% with exact 
waveforms^] so that more than 90% of the potential events can be detected. In contrast, the corresponding 2.5PN 
Taylor approximants can only detect about 50% of the potential events for massive systems (at the price of large 
biases ~ 15%). Later studies have confirmed the performance of these P-approximants {27| and assessed Q their 
need in related contexts of space based interferometers like LISA. 

C. Fourier representation of inspiral signals and validity of the stationary phase approximation (SPA) 

Independent of the choice between T- and P-approximants, another desirable approximation in data analysis for 
inspiralling compact binaries is the stationary phase approximation (SPA), which is a simple, explicit analytic approx- 
imation to the Fourier transform of the time-domain chirp (see, e.g., ps[|). In fact, most work on inspiral waveforms 
(except DIS) has used only SPA approximants to the frequency-domain chirps. In the course of our P-approximant 
work we noticed a progressive worsening of the overlap between the SPA and the "exact" Fourier transform — numer- 
ically computed by a fast Fourier transform (FFT) of the time-domain signal — (see Table II of DIS) and commented 
on these 'inaccuracies of the SPA'. In the above by SPA one means not only the problem of the formal accuracy of the 
stationary phase estimate to the Fourier transform of an analytically extended, mathematical signal but also some 
issues linked to the physics, and observability, of the real signal. In particular, in DIS, we were considering templates 
which are shut off, in the time-domain, at the last stable orbit (LSO). The present paper will also consider such 
time-truncated inspiral signals. We shall discuss this point in more detail below, but the idea is that the post-inspiral 
signal (plunge + merger) will have a frequency content very different from the inspiral one (probably pushed to much 
higher frequencies). It should, therefore make sense to try to construct filters that represent as best as possible an 
inspiral signal which lasts only up to some maximum time (time- windowing) . For such signals, DIS noted a worsening 
of the usual (frequency-windowed) SPA approximation, both as the total mass of the system increases and as the 
post-Newtonian approximation order is increased, and mentioned that this worsened performance was due to the fact 
that "such systems emit many less wave cycles in the effective detector bandwidth" centered (for initial LIGO) near 



1 This statement was proven in DIS by quantifying the convergence of the sequence of P-approximants towards some 'fiducial 
exact' waveform. In the test-particle case this waveform used the known Schwarzschild's energy function E(v) and Poisson's 
numerically computed GW flux j2fj. In the comparable mass case, it was constructed by modelling the w- dependent higher 
post-Newtonian corrections to the best known analytical results. 
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/dot = 167 Hz. In this paper /<j e t denotes the frequency at which the noise power spectrum per logarithmic bin of 
the detector is the least (or equivalently the frequency at which the detector is most sensitive to a broad-band burst). 
To avoid irrelevant, uncontrolled sources of inaccuracy, DIS used the FFT of the time- windowed chirp rather than its 
SPA to generate the frequency-domain waveform and make comparisons between the T- and P-approximants. 

The use of FFT rather than SPA in DIS makes the P-approximant computationally expensive. As will be discussed 
in detail in Sec. VI, the use of SPA or similar frequency-domain representations is far less expensive. The obvious 
need to incorporate this desirable feature makes urgent and mandatory a critical investigation of the possibility of 
marrying together the excellent performance of the P-approximants to the relative inexpensiveness of the SPA without 
a serious loss in event rate. 

Recently, some issues related to the accuracy of the SPA have been investigated. For general chirps, Chassande- 
Mottin and Flandrin |3(J have studied whether the usual conditions assumed for the validity of the SPA are necessary 
and sufficient and attempted a quantitative control of the approximation. Droz et al ]3l|] have examined other issues 
related to the accuracy of the SPA of particular relevance to gravitational wave data analysis. Unlike DIS, by SPA, 
Droz et al imply only the stationary phase estimate of the Fourier transform and discuss separately the issue of 
windowing — the fact that the signal in the time-domain lasts only from i m i n to t max or a time-window. To improve 
the SPA estimate of a Newtonian chirp, they compute the next order contribution^] (to the Fourier integral) by the 
method of steepest descent, show that it is of order relative to the leading order SPA estimate and conclude 
that it is small enough to be justifiably neglected. [Here vm is an invariantly defined 'velocity'^ related to the 
instantaneous gravitational wave frequency F and chirp mass M. by vj^ = {itMF) 1 /' i . The chirp mass is related 
to the total mass m — mi + rri2 and dimensionless mass ratio r\ — ( m ™^_^^ by M — rf^m.] They point out the 
importance of windowing, estimate the amplitude and phase modulations induced in the frequency-domain by the 
time-window and conclude that in all cases these modulations have negligible effect on overlaps. However, their 
analytic expression for the effects of windowing is only valid for values of frequencies well away from the boundaries 
of the natural frequency-window induced by the time-window, denoted by F m i n = F(t m i n ) and F max = F(t max ) - 
the gravitational wave frequencies at times t m - ln and i max , respectively. In this paper we provide a formalism allowing 
one to compute analytic approximations to the Fourier-transform of a time- windowed signal in the most crucial edge- 
frequency- domains f ~ F min and / ~ F max (including / < F m ; n and / > F max ). As first noticed in DIS and discussed 
in detail in this present work, the effect of window oscillations on overlaps (claimed to be negligible in 31 1) starts to 
be noticeable when the total mass m ^ 13AfQ and becomes very significant for m > 20 Al q. [Here we consider equal 
mass systems n = 1/4.] Since the difference between the statements in DIS and Droz et al [pll can be disconcerting 
and a serious source of confusion to the potential user community, we discuss this in further detail next. 

In DIS, what was meant in Table II (the only place where it was used) by "stationary phase approximation" was 
the product of the usual SPA by a simple Heaviside step function 9(F max — f) i.e. h(f) was truncated above a Fourier 
frequency / = -F max where F max is the instantaneous gravitational wave frequency at which the time-domain signal 
is itself terminated, assumed to be (in DIS and here) the frequency at the last stable orbit Flso- [In the following, 
we shall, for brevity, refer to this frequency Windowed Usual Stationary Phase Approximation as the 'uSPAw'.] We 
were motivated to do this from the stationary phase result itself. The SPA (to the Fourier transform of the chirp) 
says that the dominant contribution to a certain Fourier amplitude h(f) comes from a neighbourhood of time (in the 
Fourier integral) when the instantaneous frequency F(t) numerically reaches the corresponding Fourier frequency /. 
It is therefore to be expected that the signal essentially terminates at f — FLso that there is no significant power 
in the Fourier transform of the signal beyond -FLso- This is indeed true in the first approximation, as is evident from 
Fig. H below, which shows that the power in the exact Fourier transform of the time-windowed signal [computed via 
a discrete Fourier transform (DFT)] falls off much faster than the SPA for / > Flso- Moreover, as is discussed in 
detail below, in the relativistic case the usual SPA breaks down at Jlso and cannot be meaningfully extended for 
/ > Flso ■ Hence the values quoted in Table II of DIS were obtained by computing the overlap of the DFT of the 
truncated time-domain waves with the truncated SPA representation of the wave. 

On the other hand, a critical examination of Droz et al |nj reveals that their claim regarding the adequacy of the 
SPA in fact has only a restricted domain of validity. It is relevant to SPA considered as a mathematical algorithm 
to be applied to a generic smooth signal and low mass binaries (m < 13M Q ). As acknowledged by the authors, 
they do not address physical issues related to an eventual time-domain cut-off of the signal at FLso- What they 
call "Newtonian signals" are unphysical, formally defined chirps whose instantaneous frequencies are extended to 



2 We shall give below the general result for any chirp. 

3 Note that, following DIS, we shall use v = {itmF) 1 ^ , instead of vm = ^^v, in all our analysis. We also use units such 
that G = c = 1. 
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^max = -FNyquist 3> ^lsOj i n fact, better described as 'analytically extended Newtonian signals'. It is the SPA of this 
formal, analytically extended signal which is shown to produce overlaps with exact FFTs better than 0.99 even for 
massive binary systems of chirp mass Ai — 10 M , corresponding to a total mass of m ~ 23M for an equal mass 
system (77 = 1/4). These large overlaps, in our view, are not a proof of the validity of the SPA to compute, physically 
relevant, accurate frequency-domain inspiral templates, as they do not address the important issue of inspiral-signal 
termination at or near the -FLso when Flso ~ /det, the frequency at which the broad-band noise of the detector is the 
least. It turns out that for binary systems of total mass m ^ 28M© the power in the Fourier domain beyond / = Flso 
for a relativistic signal, is a significant fraction (> 10%) of the total power. If the usual (frequency-windowed) SPA 
is used in constructing frequency-domain inspiral waves we are risking the loss of more than 30 % of the events from 
binaries with masses > 28M©. [This will be illustrated in Fig. [j] below.] This is in addition to the losses induced by 
the inaccuracy of the post-Newtonian waveforms and the discreteness of the bank of templates used in data analysis. 



D. Massive black hole binaries and first detections in LIGO/ VIRGO 

Let us first establish our notation. We define the Fourier transform (FT) h(f) of a time-domain signal h(t) by 



h(t) 



dfe 



-2mft 



Kf); Kf) 



dte 27rift h(t). 



We write the (suitably transformed) output of the detector as 

Vit = h(t) + nit) , 

where h(t) is the signal and n(t) the noise. The correlation function of the noise reads 



n{t 1 )n{t 2 ) = C n (t 1 -t 2 ) 



dfS n (f)e 2mntl - t2) 



(1.1) 



(1.2) 



(1.3) 



where S n {f) = S n (—f) is the two-sided noise power spectral density. In all the present work, we shall consider a noise 
curve of the type expected for initial interferometers. For initial LIGO we take |32|, 



Sn(f)= 



So 



f 



,/>/., 



OO, f < fa 



(1.4a) 
(1.4b) 



with f s _ = 40 Hz, f = 200 Hz and So = 1 -47 x 10~ 46 Hz" 1 . In the above we have included a factor of 1/2 
[S'° ne_sldcd = 25* wo_sldcd ] because Eq. (1.4a) gives the two-sided noise; the one-sided noise would be given by 
the same formula without the factor of 1/2. The minimum of S n (f) is at / = fo and is equal to 5 m i n = 2.55o- 
However, a physically more relevant quantity is the minimum of the dimensionless quantity /i 2 (/) = fS n (f) (effective 
GW noise, see below). This is reached at the characteristic detection frequency f = f^ ct = 0.8347/o, and is equal 
to (/1™ 111 ) 2 = 2.2761/05*0. The above numerical value for f and So leads to /dot = 167 Hz and corresponds to 
h™ ia = 2.5868 x 10~ i2 . 

For VIRGO on the other hand, the corresponding noise curve is given by 



Sn(f) 



So 
2 



10 d 



f>fs, 



= 00, / < fs, 



(1.5a) 
(1.5b) 



In this case, f s = 20 Hz, f = 500 Hz while S = 3.24 x 10~ 46 Hz" 1 @. The minimum of h n (f) = y/fSjf) is 
reached at / = 103 Hz and is equal to h™ iu = 4.2902 x 10~ 22 . It should be noted that the VIRGO noise curve is 
used only in this Section, while discussing Fig. [l]. In the rest of the paper and all the Figures and Tables, the scalar 
product is defined using the LIGO noise curve. 

Anticipating on formulas to be discussed in |IIB|, the square of the signal to noise ratio (SNR) is given by 



(k, h} 2 
(k, k) 



(1.6) 
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where the scalar product is defined by 

{k ' h) = L df ^MT- (L7) 

Here, h denotes the exact signal, and k the filter used in the data analysis. We assume in this paper that the signal 
h is given by a time-truncated adiabatic inspiral signal. [For simplicity, we consider in this subsection Newtonian 
waveforms, and we approximate the Fourier transform of a time-truncated Newtonian signal by the very accurate 
improved Newtonian stationary phase approximation (inSPA) to be constructed below.] In computing p 2 we average 
over all the angles (determining both the detector and the source orientations), and we place the source at a fiducial 
distance of 100 Mpc. [Note that a coalescence rate of 10 -5 per galaxy and per year implies that in two years one 
event should happen within 100 Mpc] 

In most of the literature one uses as Fourier-domain filter k{f) the frequency-windowed usual stationary phase 
approximation (uSPAw) to estimate the SNR for an inspiral signal. W e illu strate in Fi g. [l| the loss in signal strength 
extracted by using as filter the uSPAw in LIGO and VIRGO [cf. Eqs. |L4a|) and ([Dial )], instead of using the optimal 



filter k = h (leading to the optimal SNR p 2 — (h, h)). The plot also shows on the top horizontal axis the last stable 
orbit frequency corresponding to the total mass in question. The left vertical axis shows the SNR extracted and the 
right vertical axis shows the sensitivity, h^ 1 = [/S„(/)] _1 , (both of which are dimensionless) of LIGO-I and VIRGO 
instruments. While reading the sensitivity curve one should use the top and right axes and while reading the SNR 
curve one should use the bottom and left axes. The SNR values plotted in Fig. [l] have been computed numerically 
by inserting the relevant values of h(f) and k{f) in Eq. ([O]). Though we did not use it, it might help the reader to 
see the analytical expression of p 2 obtained in the simple approximation where k(f) ~ h(f) ~ h uspa,w (f). Using the 



equations of Sec. O below (and averaging over angles as explained in Sec. IV A below) leads to 



{p2) * ifc i7J I 7^)7^(7)' (L8) 

where (..) den otes the angular average, d the distance to the source, v(f) = (irmf) 1 ^ 3 and S n (f) the two-sided noise 
given by Eq. ( 1.4a ). We indicated no precise detection threshold in Fig. |l| because this depends on many parameters 



(like the number of detectors involved). The reader should, however, have in mind that a reasonable detection 
threshold is, at least, Pthreshoid ~ 5. 

We note that the effective sensitivity of LIGO-I peaks near a frequency of 167 Hz which is the last stable orbit 
frequency for a binary of total mass of about 27M Q . The effective sensitivity of VIRGO peaks at a much lower 
frequency of 103 Hz. This low-frequency sensitivity of VIRGO means two important things: Firstly, lighter binaries 
(i.e., m < 30Mq) are integrated for a longer time in the low frequency regime and, therefore, the corrections to 
the Fourier transform introduced in this paper are less important for such systems. This means that the uSPAw is 
quite good in extracting the full signal power of such binaries as evidenced in Fig. [|. Notice that, for VIRGO, the 
uSPAw curve follows the inSPA curve for m < 30M Q . On the contrary, LIGO's lower sensitivity to lower frequencies 
makes it important to include the corrections to the Fourier transform of LSO-truncated signals from binaries of mass 
m > 15M Q . In LIGO's case uSPAw extracts only 75% of the full SNR implying a loss of more than 40 % of all massive 
binary coalescences. On the other hand, the low-frequency peak of VIRGO sensitivity means that we will have to 
employ the accurate Fourier domain models discussed in this paper for more massive binaries, i.e. m > 30Mq. [Note, 
however, that the low-frequency sensitivity of VIRGO means that, for low-mass and medium-mass binaries, it is even 
more crucial to use P-approximants (instead of the usually considered T-approximants) than for LIGO, in order to 
accurately keep track of the phasing of the many cycles accumulated at low frequencies.] 

It is fair to say that at present the most well-understood gravitational waveform is the inspiral one and thus the 
only reliable templates correspond to inspiral signals. It is also generally believed that binary black holes are better 
candidates for gravitational wave sources than binary neutron stars due to their larger masses (the average mass of 
observed black hole candidates is around 8M Q [M). Theoretical computations based on stellar evolution indicate 
that binary black holes with individual masses < 15M Q may be the only known sources that exist (hopefully) in 
sufficient numbers |35| -^8|] . When looking at Fig. [y, one clearly sees the importance of dealing with binary black holes 
with total masses in the range of 28 — 30M©. They lead to signals with the best SNR. However, it is precisely for such 
systems that the FLso is around the middle of the detection bandwidth for initial LIGO i.e. -FLso ~ /det- The most 
likely sources to detect are in the problematic region discussed in this paper. This makes it imperative to not lose 
SNR when dealing with such signals and provides the other major motivation for our work. If our assumption that 
the best models of inspiral waveforms must be abruptly shut off in the time domain holds, it is essential to use the 
improved SPA formulas discussed in this work in order to maximise our chances of detecting inpiralling binaries. The 
analysis presented in this paper provides insights and techniques to deal with binary black hole signals in probably 
the most crucial mass range. 
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E. Summary of the present paper and proposals for data analysis groups 



In this paper we propose analytical approximations to the Fourier-transform of the LSO-truncated time-domain 
inspiral waves that are very accurate even for the massive black hole binaries, the most likely sources for LIGO and 
VIRGO [overlaps with FFT > 0.99 for to < 40M Q ] and are at the same time computationally inexpensive. We call 
our final new, frequency-domain filters the SPP approximants because they combine the computational convenience 
of stationary phase approximants with the accuracy of the (time-domain) P-approximants. Our strategy is two-fold: 
On the one hand we introduce a correction factor C(f) to the usual SPA for / < / up < -Flso which improves the SPA 
by taking into account the "edge" oscillations present when / < / up . (/ up will be defined below. For Newtonian-like 
signals / up = -FlsO) while for relativistic signals / up < -Flso-) On the other hand, we introduce a new approximation 
to the Fourier transform for / > / up which efficiently recovers the signal power around and beyond the frequency 
corresponding to the last stable orbit. These features are important new steps forward as there was no formalism 
until now that could compute (especially for relativistic signals) Fourier transforms analytically for / ~ -Flso > an d in 
particular / > -Flso- These new features now make it possible to generate templates directly in the Fourier- domain, 
leading to a saving on the computational cost of template generation by a factor of 10 or more. 

Our concrete proposal to the interferometer data analysis groups that are building the gravitational wave search 
software and wish to have Fourier-domain filters which are both accurate and fast-computed, is thus the following: 
First, we confirm that for accurate post-Newtonian template generation of binary systems of total mass to < 40Mq 
one needs to use a frequency-domain version of the P-approximant (previously defined only in the time-domain) . For 
to < 5M Q a straightforward (uncorrected for edge-effects) SPA of the P-approximants is sufficient. (They match 
with the exact DFT of the same time-signal with overlaps > 0.999.) On the other hand, in the total mass range 
5Mq < to < 40M Q , and assuming that one wishes an accurate frequency-domain (f-domain) representation of a 
time-windowed signal it is crucial to use our new SPP approximants. For m > 40M© a straightforward DFT is 
recommended (but, anyway, the signal is not known with enough precision in this high mass range, where the plunge 
and merger signals become observationally important). 

It is important to stress the position we assume in this paper: Given the absence of any detailed and precise 
information about the plunge signal today, we suggest that a time-truncated chirp (time-windowed signal) is currently 
our best bet and the modified SPA presented in this paper is the appropriate Fourier-domain representation one must 
use. However, this should not be taken to imply that we are claiming to have logically excluded the other possibility 
that the f-window may turn out to be the better choice, when we have further details about the transition from inspiral 
to plunge and, about the plunge waveform. Even so, we emphasize that a definitive contribution of the present work 
is to provide explicitly for the first time the frequency-domain version of the time-domain P-approximants which were 
shown in DIS to bring indispcnsible improvements over the usually considered T-approximants. Consequently, even 
in the unlikely case where a straightforward frequency-window turns out to be a better model than the time-window 
assumed in most of this work, one will still require the formulas given in this paper [with the trivial change of replacing 
the correction factors C(£) by a 9 function #(-Flso — /)] to generate sufficiently accurate f-domain filters. Thus this 
work may not be the complete final answer but only a step ahead and a partial contribution towards defining good 
f-domain filters. Assumptions that seem the best we can accept, require special tools for their analysis and this paper 
provides them. 

This paper is organised as follows: In Sec. II A , II B , [I C as a prelude to later technical material, we introduce several 
useful physical notions and emplo y the m to give a preliminary discussion of the questions raised by the detectability 
of massive-binary signals. In Sec. |ll D| wc summarise the mathematical tools used in the p aper to estimate the time- 
truncated chirps. In Sec. pi we c onsider time-windowed Newtonian- like signals. Sec. [II A provides a short summary 
of th e well known SPA. Sec7|lIB sets up the basic equations to discuss the Fourier transform of time- windowed signals. 
Sec. Ill C estim ates th e e dge co ntribution to the Fourier transform coming from the non-resonant integral. This is 



followed by Sec. Ill D and HIE where we elaborate in detail the construction of optimal analytic approximations to 



the Fourier transform of the time-windowed gravitational wave chirp (improved Newtonian SPA). I n Sec . Ill F we 
compare and contrast in detail the usual SPA with our improved SPA for Newtonian-like sign als. S ec. IV A addresses 
the new issues related to the Fourier-transform of time-windowed relativistic signals. In S ec pVB] we present a new 
method to estimate the small non-resonant contribution in the relativistic case. In Sec. IV Q we construct a new 
form of improved SPA for such signals (improved relativistic SPA). Combining this improved relativistic SPA with 
the P-approximants of DIS leads to the construction of the frequency-domain SPP approximants. In Sec. [v] we use 
the SPP approximants constructed earlier and investigate their faithfulness and effectualness in detail for the test 
mass case. Based on this we comment on the corresponding situation in the comparable mass case. In Sec. wc 
compare the computing costs for template generation using the time-domain FFT with corr esponding costs for the 
frequency-domain SPA and improved SPA both for Newtonian and relativistic cases. Sec. VII contains our concluding 
remarks. 
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II. PRELIMINARY DISCUSSION 



As a preface to the technical treatments of the following Sections in which we shall construct optimal analytic 
approximations to the Fourier transform of the gravitational wave inspiral signal h(t) let us start by discussing some 
general issues which are central to this paper. 

A. Wiener niters and time-truncated inspiral signals 

We briefly recall the principle underlying the optim al l inear filter technique (Wiener filter) . A (real) linear filter is 



a linear functional of the detector's output, h au t, Eq.(1.2), say 



K[h out }= dtK(t)h aat (t)= dfK(-f)h out (f)= dfK*(f)h out (f). (2.1) 

J — oo J — oo J — oo 

Let us associate to any K(t) the time-domain function k(t) such that its FT equals k(f) = S n (f)K(f) and let us 
introduce the Wiener scalar product (defined on real time-domain functions) 

/OO J J* pOC pOQ 

Sjfj9*WHf)= I I dt 1 dt 2 g(t 1 )wi(t 1 -t2)h(t 2 ), (2.2) 



oo J — OO 



where, w x {r) = J°° ^ e 2 "^ = ^(-r), (2.3) 

is the convolution inverse of the noise correlation function C„(t) = C n (—r) i.e. 

( Wl *C n )(t)=5{t). (2.4) 

[Here * denotes the convolution product: (g * h)(t) = dr g(r)h(t — r)]. With this notation the action of the filter 
K on /i out reads 

K[h out ] = (k,h out )=S + N, (2.5) 
where S is the filtered 'signal' and TV the filtered 'noise', defined by 

S = K[h] = (k,h); N = K[n] = (k, n) . (2.6) 



The definition of the symmetric Wiener scalar product Eq. (2.2) is such that the statistical average of a product of 



filtered noises simplify: (k\ , n) (k 2 , n) = (ki,n)(n,k 2 ) = (feijfea). In particular, the variance of the filtered noise N 
reads N 2 — (k, n) 2 = (fc, k), so that the square of the signal-to-noise ratio (SNR) for the filter defined by any function 
k{t) reads 

where we have defined the 'overlap', or normalized ambiguity function, between k and h 

Q{k.h) = —= . (2.8) 

y/(k,k)(h,h) 

Schwarz's inequality guarantees that \0(k, h)\ < 1, the equality being reached only when k(t) = Xh(t). (We work here 
in the s pace of real signals.) For a given signal h(t), the choice of filter K <-> k which maximizes the SNR p is, in view 



of Eq. (2/7), k{t) = \h(t), where the proportionality constant is unimportant and can be taken to be one ('Wiener 
Theorem'). This optimal linear filter theorem applies when the full time-development of the signal h(t) is known, the 
noise is stationary and has a known spectral distribution S n (f). 

It is important to note that, saying that the best 'associated' filter k(t) is simply k(t) — h(t), means that the best 
time-domain filter K(t), which must be directly correlated with the detector's output, iT[/i ou t] = dt K(t)h ou ±(t), 
is a non-local (in time) functional of h(t). Explicitly, 
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K(t) = (wi *h){t) = / dTWi(r)h(t-T). (2.9) 



This poses the question whether, for an off-line analysis of the data, one would like to store a bank of these non-local 
time-domain filters K(t), which densely cover the expected parameter space. In the present paper, we shall assume 
that one anyway computes (nearly) on-line the Fast Fourier Transform (FFT) of the detector's output (which is 
anyway needed to factor out the frequency dependent effect of the interferometer on the GW signal), and we shall set 
ourselves the task of providing the best possible analytical representations of the Fourier transform of the expected 
signals h(f). Moreover, the availability of a fast Fourier algorithm makes the filtering problem computationally less 
intensive in the Fourier-domain. It is well-known that the computation of a di screte correlation for all (discrete) 



time lags between the output and the filter [ie the discrete version of Eq. (2A) requires 0[N ] operations in the 



time-domain while it takes only Q\N log?( iV)1 operations in the Fourier domain (because a time lag r adds a factor 
cxp(27tz/t) in the f-domain version of Eq. ( |2.l| ), which is equivalent to computing a certain inverse Fourier transform). 
Discrete correlation in the f-domain suffers from spurious correlations for non-zero time lags but this is easily taken 
care of by padding the tail part of a template with large number of zeroes (see, eg. Ref || for details.) 

As the problem of the locality/non- locality in time will be crucial to our discussion of the inspiral signals, we wish 
to give an alternative discussion of Wiener's optimal linear filter theorem. Indeed, another proof of the theorem can 
be obtained by introducing a 'whitening' transformation say mi , which simplifies the properties of the noise. We 
define the 'whitening- kernel' aii(r) by 

wi (r) = / _L=e 2 ^ ; Wi (/) = -== . (2.10) 



For any function g(t), the action of the kernel wi on g (i.e. the 'whitened' version of the function g) can be denoted 
by 



gi(t) = (wi*g)(t)= dTwi(T)g(t-T). (2.11) 

J — OO 

The name 'whitening kernel' conies from the fact that the transformed noise ni(t) is 'white' i.e. uncorrelated: 

roc 

Tu(ti)ni(t 2 )= / dfe 2 * if ^-^ =5{t 1 -t 2 ). (2.12) 



Note that wi is simply the convolution square-root of the Wiener kernel w\ introduced above: wi * wi = wi. 
The Wiener theorem states then that, after having whitened all the functions, the optimal filter is simply the usual 
straightforward correlation between the (whitened) output and the (whitened) signal i. e. 

/oo 
dth$(t)hT(t), (2.13) 

where hi = wi * h, h° ut — wi * h out . In other words, we can think of the optimal filter as being local-in-time 

2 2 2 2 

after the application, to all signals, of the convolution kernel wi. When working with the transformed time-domain 
functions hi(t), JW ut (i) • • •, we shall say that we work in the 'whitened time-domain'. In this language, the best filter 

2 2 

in the whitened time-domain is to simply correlate (as when trying to visually superpose two time-functions) the 
output with a copy of the signal. This 'whitened time domain' is conceptually useful in the present context because 
it introduces only a small non-locality (by small we mean here much smaller than the non-locality introduced by 
the Fourier transformation). Indeed, as the function l/y/S n (f) has a rather flat maximum, its Fourier transform 
wi(t) is nearly a delta function as seen in Fig. || wherein we have plotted the whitening kernel for the initial LIGO 
interferometer. More precisely, wi(t) is an even function made of a positive spike around r = 0, followed (on each 
side) by a slightly negative wing which decays fast towards zero as |t| — > +oo. The half- width at half-maximum of 
the central spike is 0.18 ms. The location of the wings is around t = ±0.002 s. Therefore the non-locality contained 
in the whitening transformation is only between 0.2 ms and 2 ms (depending on the function on which it acts). 
This non-locality together with the one and a half cycle in wi is sufficient to efficiently damp both high and low 
frequencies so that hi — wi * h is a chirp whose amplitude is important only when the instantaneous frequency is 
around f p — 165 Hz (see below). We shall use below this whitened time domain picture to discuss the important 
features of the expected chirp that we should try to model as well as possible. 
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In the present paper, we shall be primarily interested in massive compact binaries with total mass m = mi + TO2 
in the approximate range 3M© < m < 4OAf . We recall that the GW signal from a compact binary is made of an 
inspiral signal, followed, after the last stable orbit is reached, by a plu nge sig nal which leads to a final merger signal. 
Thanks to the analytical work on the motion || of and GW emission [pi| , [39[ from general relativistic binary systems 
we have quite a good analytical control of the inspiral signal. In the present paper, we shall further argue that we have 
also a rather good analytical control on the location of the last stable orbit (LSO), i.e. on the transition between the 
inspiral and plunge. First, DIS [ fl3[ introduced a new, more robust approach to the determination of LSO based on 
the invariant function e{v). More recently |4C|| , a new approach to the dynamics of binary systems has confirmed the 
result of DIS (saying that the LSO was slightly more 'inwards' for comparable mass systems than as predicted by the 
test-mass limit) and predicts values for the important physical quantities at the LSO (notably the orbital frequency) 
which are even nearer to the (Schwarzschild-like) ones obtained in the test-mass limit. We anticipate that further 
analytical progress in the problem of motion of binary systems ( fli| , combined with the LSO-determination techniques 
given in p3| and pOp^ ] will soon allow one to know with more certainty and more precision the gravitational wave 
frequency at the LSO. Pending such a determination, we shall use as fiducial value for the GW frequency at the LSO 
— when we need it for simple analytical estimates — the usual 'Schwarzschild-like' approximation 



Flso = 4400 



Me 
m 



Hz. 



(2.14) 



However, in our actual numerical calculations and plots we shall use the ^-dependent -FLso corresponding to the 
approximation used for the energy function Ea(v). For instance, in the case of the 2PN P-approximant P4 we have 



^LSO 



4397.2 



2 - 



1 



16 



L v 2 



:!6 



Mr. 



Hz. 



(2.15) 



In the equal-mass case (77 = 1/4) this yields Flso = 5719.4M /to Hz. Note that the most recent dete rmin ation of 
the LSO @ suggests that when 77 ^ the GW frequency at the LSO lies between Eq. (|2.14| ) and Eq. (|2.15|) : 



Flso - 4397.2(1 + 0.31557?) { ^ 

TO 



Hz. 



(2.16) 



Preliminary studies indicate that the plunge signal, emitted during the fast fall of the two masses towards each 
other following the crossing of the LSO, will last (when Ar) ~ 1) only for a fraction of an orbital period (see an d 
Hg. As usual, one can also assume that the subsequent merger signal linked to the formation of a black hole of total 
mass ~ m contains significantly higher frequencies than the inspiral ones. Indeed, the characteristic frequency of the 
merger signal may be taken to be given by the real part of the most slowly damped quadrupolar normal mode of a 
black hole (which when neglecting the black hole spin, has a complex circular frequency wiTObh = 0.37367 — 0.08896 i 
PU), i.e. (with TObh ~ m) 



fn 



bh 



0.374 

27TTO 



12000 — Hz . 

TO 



(2.17) 



Eqs. ( p. 14 ) and (2.17) lead us to accept that there is a significant frequency separation between the inspiral and plunge 
signals and the merger one: /merger/^LSO ~ 2.75. Therefore, if we restrict our attention to systems such that the 
characteristic detection frequency /d e t defined by the noise curve, stays logarithmically nearer to Flso than to /merger, 
it seems plausible that a good filter to use for GW detection can neglect the (ill-known) merger signal, but shoul d try 
to model as accurately as possible the inspiral and plunge signal. For the initial LIGO noise curve, Eq. ( 1.4a ), the 
characteristic detection frequency /det is 167 Hz. It is then for a total mass to < 43.5M Q that /det/^Lso < /bh//det- 
We shall see below that, just before reaching the LSO, the inspiral signal is still significantly 'quasi-periodic' (with 
> 6 cycles before a significant change in instantaneous frequency). By contrast, though the plunge signal may not 
decay monotonically and may be oscillating, it seems reasonable to assume that the plunge lasts only for a fraction of 



4 We assume here that we are in the generic case where the spins of the coalescing objects are smallish compared to their 
maximal Kerr value. 
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the orbital period Tlso = ^-F^go- Thus, in the absence of a precise knowledge of the plunge signal, a good model of 
the time-domain signal consists in abruptly shutting off, by a step function 6>(£lso — t) the (adiabatic) inspiral signal 
beyond the time t^so when the last stable orbit is reached. We also tested formally the robustness of our approach 
by showing that our model above has a good overlap with a signal which decays smoothly on a time scale of a few (up 
to 3) -Flso beyond the LSO. Because of the likely oscillatory behaviour of the plunge signal, details of the oscillations 
are necessary for any further improvements and we are currently working towards improving our understanding of 
the transition between the inspiral and the plunge p3[] . 

These considerations motivate us to propose that, in the absence of knowledge of the optimal filter which should be 
fcoptiroal(^) = k cxac t(t), our best bet is to use as (sub-optimal) filter the time-truncated inspiral signal /iinspirai(i)^(^LSO — 
t). In other words, we think that the best strategy is to use all the information available, in the time-domain, about 
the signal and to replace the transient plunge and higher frequency merger signals by zero, as a measure of our current 
ignorance. But having settled on this tactic in the time-domain, the aim of this paper is to provide the best possible 
frequency-domain description of such a time-windowed signal. We shall see in detail below that the Fourier transform 
htvr(f) = FT[/itw(i)] of the time- windowed signal ht w (t) = /iinspirai(0^(^LSO — t) has a non-trivial structure which is 
not captured by the usually considered frequency-windowed stationary phase approximation h spllw (f). In particular, 
for massive systems, a significant fraction of the total power is contained in the 'tail' of h tw (f) beyond / = Flso- 
The general result Eq. ( |2.7| ) for SNR obtained with any filter then says that the Fourier-domain filter /itw(/)j though 
sub-optimal, should be still significantly better than /i sp aw(/) because (under our assumptions about the plunge + 
merger signal) it has better overlaps with the exact signal. We shall return below to this important issue and give 
further arguments (in the time-domain) to confirm the superiority of h tw (f) over 

fr S paw(/) (see Fig. | and Fig. [h] and 

text around it). 



B. The number of useful cycles 



Often one mentions that, in the total time development of an inspiral signal, the total number of gravitational wave 
cycles 

N tot = -(4 nd - ^ Cgin ) = dF [- , (2.18) 

is very large. Here 4> is the gravitational wave phase, </> on d is the phase at the end of the inspiral regime (defined 
by the last stable orbit for sufficiently massive systems, i.e. for the black- hole-neutron-star and the black- hole-black- 
hole systems), while 0bcgin is the phase when the signal enters the lower frequency (seismic) cutoff of the detector 
bandwidth. We have also rewritten N tot as an in tegral over the running instantaneous gravitational wave frequency 
F. However, the large number A^ot, Eq. ( |2.18| ), is not significant because the only really useful cycles are those 
which contribute most to the signal-to- noise ratio (SNR). To have a clearer idea of what one might want to mean 
by the notion of a useful number of cycles, let us first introduce the instantaneous number of cycles spent near some 



instantaneous frequency F. It is naturally defined by multiplying the integrand in Eq. (2.18) by F, considered as the 
length of an interval ±AF = ±-j around F i.e. 

N < F ^i§ s 7^Tf < 2 - 19 > 

where we have used d(f>/dt = 2nF(t). Note that the instantaneous N can be considered either as a function of the 
running frequency F(t), or, directly, of time. 

The instantaneous number of cycles plays an important role both in defining the observability of a signal, and in 
controlling partially the validity of the stationary phase approximation. The square of the optimal SNR reads 

s\ 2 _ f + °° jr IM/)I 2 



* 2s ^j -y_ *m- ,2 - 2o) 

In the stationary phase approximation (discussed at length and improved below; but here we use standard results for 
orientation) the modulus of the Fourier transform of the real signal h(t) = 2a{t) cos </>(£), reads |/i(/)| — a(tf)/ \f F(tf) 
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where tf is the time when the instantaneous frequency F(t) reaches the value /^] Therefore, the squared modulus 
can be written as 



« 2 (/) _ 1 

V l \J )\ - 

Finally, the SNR can be rewritten as 



^-mi^p N{f)a2{f) - (2 - 21) 



2 = flV - f +0 ° # N (f) a2 (f) [ + °° d ± MiZl (9 99, 

p -{n) -J^ f fs n {f) ~J-oo f hUf)' [ - ] 

where we have introduced the notation h 2 s {f) = N(f)a 2 (f) and h^f) = fS n (f). Here, h^{f) is the usual squared 
amplitude of the effective gravitational wave noise at the frequenc y /, i.e. the minimum gravitational wave ampli- 
tude detectable in a bandwidth ±//2 around frequency /. Eq. ( |2.22 ) exhibits that the squared amplitude of the 



corresponding effective gravitational wave signal is h 2 (f) = N(f)a 2 (f) = f 2 \h(f)\ 2 , i.e. that the "bare" amplitude 
a(f) ee a(t(f)) is effectively multiplied by y/N(f) §j|. Eq. ( |2.22| ) also exhibits the relative weight with which each 



cycle counts for detectability purposes. Per logarithmic frequency interval this weight is simply 

w{f)^a 2 {f)/h 2 n {f). (2.23) 
Therefore it is natural to define the number of useful cycles as 

N useiul ee ( I*'™ % w(f) N(f)) ( - ,,, / ) I . ( 2.24 ! 



'.2 




where F min is the low- frequency seismic-cutoff below which ft.„(/) is essentially infinite and where the upper-cutoff 
-Fmax is the frequency at which the signal itself shuts off. For illustration, we list in Table || the number of useful 
cycles and the total number of cycles for representative systems and orders of approximation. For Newtonian chirps 
the total number of cycles is 



. Nl . Vll (7rm/ B ) 5/3 - (7rm/ max ) 5 / 3 



NZt m = 32 ^ ' (2 ' 25) 

where f s is the seismic cutoff. The total number of cycles for relativistic chirps is alw ays sm aller than wt . From 
Table | it is clear that the number 2V usc f u i is usually much smaller than N tot , Eq. ( 2.1 8| ) . Note that for massive 



systems iV use f u i becomes quit e small. The number of useful cycles given in Table | have been c ompu ted for the initial 



LIGO noise curve Eq. ( 1.4a ). The corresponding numbers for the VIRGO noise curve Eq. ( 1.5a ) would be larger 
both because the VIRGO sensitivity curve peaks at a lower frequency, and because it is flatter. 

To explore in more detail the case of systems with a small number of useful cycl es we display in Fig. ^| (c 



linear-log plot) the various factors of the logarithmic integrand on the RHS of Eq. ( |2.22 ) for two different binary 



systems. The instantaneous number of cycles N(f) in the Newtonian approximation is plotted, together with square 
of the amplitude a 2 (f), their product the effective gravitational wave amplitude h' 2 (f) = N(f)a 2 (f), the reciprocal 
of effective noise h 2 (f) — fS n (f), (cut off after F max = Flso) an d the power per log bin of the square of the SNR 
dp 2 /d(log /). On the left, the scale on the y-axis corresponds to N(f). On the right it corresponds to the amplitude on 
an arbitrary scale. Other quantities are on an arbitrary scale. The top panel is for a lighter mass binary (mi = IAMq, 
n%2 = 10M Q ) and the bottom panel for a heavier one (mi = m 2 = 1OM ). 

The last Figure exhibits two useful lessons which are well-known but are particularly important to keep in mind 
when reading the present paper. First, because of the mass-scaling of the gravitational wave frequency at the last 



stable orbit, given in the lowest (Schwarzschild-like) approximation by Eq. (2.14), it is only for systems with total 



mass m = mi + mi > 13Mq that the peak of the SNR logarithmic frequency- distribution f p becomes comparable, 



5 In the following, it will be necessary to distinguish carefully between the instantaneous frequency F and the Fourier variable 
/. In the present Section this distinction is not very important and we shall freely change notation / <-» F. Similarly, the 
gravitational wave flux and factored flux functions which following standard notation was denoted by F(v) and /(«) in jl^] 
are denoted by J-(v) and l(v) in this paper to avoid confusion with the instantaneous gravitational wave frequency and Fourier 
variable respectively. 
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within a factor of two, to F max — -Flso- This statement critically depends on the characteristic fre quen cy entering 
the considered noise-curve. For instance, in Fig. || we have used the initial LIGO curve Eq. ( 1.4a ) for which 



fp = 0.825/o = 165 Hz. Note that the peak of the logarithmic SNR integrand, f p , is very close to the minimum 
of the effective noise amplitude h^(f) = fS n (f), which is located (as mentioned above) at /dct = 0.8347/o — 167 
Hz. This is because, the frequency dependence of the factor N(f) oc (nv 5 )" 1 oc /~ 5 / 3 (which in the effective 
signal h 2 s (f) = N(f)a 2 (f), favours lower frequencies) is nearly compensated by the frequency dependence of the bare 
amplitude a 2 (/) oc v 4 oc / 4 / 3 (which favours higher frequencies). 

A second lesson, to be drawn from Fig. || is that the number of useful cycles also becomes small in the same 
problematic case of massive systems. To see this more clearly, let us write down the explicit expression for the 



instantaneous number of cycles. In the Newtonian case (for which the basic formulas are recalled in Sec. Ill A below), 
one has 

5 11 

-^Newtonian (/) = 7T^— "j * > (2.26) 

where v — (ttto/) 1 / 3 . The lowest value of N is physically that formally reached at the upper frequency-cutoff 
/ = -Fmax = -Flso- For wlso = 1/V6 (the "Schwarzschild" value), the above equation reads 

N (F ) - 5 fiV 2 1 ~ 5 - 8477 (O 97) 

^Newtonian (, -FLSO J — O T~ — : , \ Z - Z <) 

247T 4?7 An 

where we recall that 77 < 1/4 and that the upper value 77 ma x = 1/4 is reached for equal mass systems m\ = 7712. 
Therefore for comparable-mass, massive systems, the Newtonian approximation suggests that the useful number of 
cycles will be rather small (~ 6) and concentrated near the LSO. As we shall see later, if we were interested in 
estimating the Fourier transform (FT) of an analytic Newtonian-like signal, even such a smallish number of cycles 
(and even a smaller one, down to N ~ 1) would be enough to ensure that the leading correction to the stationary 
phase approximation is small. However, the complication comes from the combination of two facts: (i) the signal 
essentially terminates at the LSO crossing time thso, and (ii) one crucially needs a relativistic description of the 



evolution near the LSO. Using the formulas and the notation of Sec. IV below, we find that the relativistic prediction 



for the instantaneous number of cycles (in the adiabatic inspiral approximation) is 

AW(/) = -^« 4 f^. (2.28) 
J7T -ryv) 

By definition of the LSO (see e.g. the discussion in DIS) the derivative E'(v) vanishes at v = vlso- Therefore, the 
instantaneous number of cycles is smaller in the relativistic case than in the Newtonian one and actually tends to zero 



near the LSO. We shall tackle in Sec. IV below the problem that this vanishing of A(Flso) causes for the stationary 
phase approximation. In this introductory Section let us only illustrate the problem by plotting the Newtonian and 
relativistic values of N(F(t)). In Fig. [| we plot the instantaneous number of cycles for the Newtonian and second 
post-Newtonian P- approximant waveforms in the last few cycles of the binary inspiral for a (20AF Q , 20M Q ) system. 
We also show the development of the waveform in this interval on an arbitrary scale. These plots demonstrate how 
the number of useful cycles diminishes as one gets close to the LSO and lead us to anticipate the subtleties in the 
detectability of signals whose LSO is near the most sensitive part of the frequency response of the detector. 



C. Loss of SNR due to edge effects 

As we already mentioned, in addition to the problem of a vanishing instantaneous number of cycles near the LSO 
in the relativistic case, the main problem with the acccuracy of the stationary phase approximation comes from the 
fact that the Fourier transform of a time- windowed signal h tw (f) = FT[/i tw (^)] differs from the frequency- windowed 
SPA h specw (f) — 6(F max — f)h spa (f) because of some 'edge effects' in the frequency-domain, linked to the abrupt 
termination of the signal in the time-domain. These edge-effects comprise some additional oscillatory behaviour in 
h(f) for / < Flso j as well as a decaying oscillatory tail in the usually disregarded frequency-domain for / > Flso- 

Let us here anticipate on the results below and use a first-order approximation to discuss the main features of the 
corrections brought by the time-windowing. Roughly (see below) the exact Fourier transform can be written as 

M/)^X(/) + £(/), (2-29) 

where h^^(f) = /i spa (/)6'(F ;l max — /) is the usually considered frequency-windowed SPA. The difference £•(/) is approx- 
imately of the form 
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£(/) ~ V(f)h B ^(f) 



(2.30) 



where, V{f) = C(f) - 6(F n 



I) 



(2.31) 



with the correction factor C(f) given by Eq. ( 3.23 ) below (with any choice of C(/) i n the present approximation) and 
where h spa (f) is some smooth continuation of h s ^(f) from the domain / < -FLso to the domain / > i'LSO- (For the 
present purpose one can think that h spa (f) is simply given by the Newtonian approximation). 



Starting from Eq. ( 2.29 ) one can compute the overlap between the exact h(f) and the usually considered frequency- 
windowed SPA h spa (f) 



o = 



h spa ) 

win ' win/ 



(2.32) 



As recalled above [see Eq. 
order in e the overlap Eq. 



^{h x ,h x ){hT 

( |2.7])1 this overlap, if it is significantly smaller than one, represents a loss in SNR. To lowest 
( |2.32| ) differs from 1 by 

1 



l-O 



2\\h 



x 



\e\\ 2 -\(E,hx)\ 2 ).. 



(2.33) 



where ||e|| 2 = (e , e), and hx = hx/\\hx\\- In inserting the explicit result Eqs.( [2.3(i| ) and (2.31) for e one sees that the 
second term on the RHS of Eq. (2.33) is much smaller than the first (because the oscillations in T>(f) are integrated 
against the smooth variation of h spa (f)). Finally, if we define the weight function 

/|/r^(/)| 2 ^ N(f)a 2 (f) hl{f) 



Sn(f) 



KU) ~ Kif) ' 



which is the full logarithmic weight function appearing in the squared SNR, Eq. (2.22), we can write 

I /' x ,//• . / /•'•••• - ,//' \ 

l-O 



(2.34) 



(2.35) 



/ w " w " \Jo f 

As will be discussed later (see Fig. ||) the function 1>(f) — C(f) — 9(F max — /) is concentrated in an interval of order 
J F(t max ) around / ~ F max and decays on both sides [like l/C(f) oc l/(f — F max )] when / gets away from i^max- 
The total integral of \T>(f)\ 2 is finite and of order unity. Thus, we see from Eq. ( 2.35| ) that when the characteristic 
frequency / p around which er(/) is concentrated satisfies / p <C F max we shall have a rough estimate 



l-O 



g(Fmax) V- F (^max) _ a(F max ) 
^(/p) -?max 



while in the opposite limit where f p 3> F max , we get the rough estimate 



l-O 



F(t n 



VN(F max ) 



(2.36) 



(2.37) 



In the case where / p -C F m ax the factor cr(F max ) / 'cr(f p ) in the RHS of Eq. ( 2.36 ) is very small. Therefore, even 
if the number of cycles iV(F max ) is not very large, Eq. (2.36) predicts that a simple frequency-windowed SPA 
will have excellent overlap with the exact h(f). On the other hand, Eq. (2.37) shows that in the reverse limit 
f p -Fmax which means in fact when / p > -F max , the overlap will become bad if yj N(F mSLX ) is not very large. As 
we have seen that \/N(F max ) gets as low as y/5.85 ~ 2.4 in the Newtonian case, and reaches smaller values in the 
relativistic case, we expect that the cases where the frequency-windowed SPA has a bad overlap with h(f ) are those 
where / p becomes comparable, say within a factor of two, to F u 



= Fi 



LSO- 



We recover the same conclusion as 
above, which was the conclusion already pointed out in DIS: namely the signal from massive systems [m ^ 13M if 
f p = 165 Hz, corresponding to /q = 200 Hz, more generally, m > 13 (165Hz// p ) M Q ], when treated (as they should 
be) relativistically will be badly represented by the usual frequency-windowed SPA. This conclusion obtained from 
an analytical approximation is borne out by the numerical computations shown in Fig. [j] above (see also Table || 
below). It is mainly for such systems that the work presented in the following Sections will be mandatory. But even 
for lower mass systems, we shall construct here for the first time the Fourier-domain version of the (time-domain) 
P-approximants introduced in DIS. Since P-approximants provide better templates than the usually considered T- 
approximants |L3], the work of this paper will be useful for all types of systems, even the less massive ones. 
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D. Fourier transform of time-truncated chirps 



To introduce the detailed analysis that we shall give in the following Sections, let us start by delineating some general 
mathematical facts about the integrals we have to deal with. We will be interested in evaluating the Fourier-transform 
h(f) of a time-truncated chirp h(t) — 2a(t) cos 4>(t)9(t max —t). After decomposing the cosine into complex exponentials, 

the Fourier integral leads to a sum of two integrals of the form f^™* dta(t) e^s ^ with phases ipf = 27r/i ± <fi(t). 
Let us then, for generality, discuss the properties of integrals of the type 



/(e) = J dta{t)e l — . (2.38) 

Here, we have introduced a formal 'small parameter' e (set to unity at the end of the calculation) to formalize the 
fast variation of the phase compared to that of the amplitude. 

Let us first note that: (i) if the phase has no stationary point tp(t) ^ for t 6 [t a ,tb], (ii) if the amplitude vanishes 
smoothly at the edge points t a and if,, which might be pushed to ±oo and (iii) if the functions a(t) and ip(t) are 
smooth (C°°) within the interval [t a ,tb], the integral 1(e) tends to zero with e, faster than any power. This can be 
seen by integrating by parts. To simplify the calculation we can [thanks to the assumption (i)] use ip as integration 
variable. This yields 

b drp A(ip)e^ , (2.39) 

where A(ip) — a(t)/ip(t) where t(ip) denotes the (unique) solution in t of tp = ip(t) and where ip a = ''Pita), ipb = 
L J t(v») 

4>(tb). Using e l ^/ £ — |^ (e l ^/ £ ), integrating by parts, and using [thanks to assumption (ii)] the vanishing of A(ip) 
at the edges, leads to 

/(e) = ie dt/jA'(ip)e— . (2.40) 

Jlpa 



Using [assumption (hi)] the vanishing of all the derivatives of A(ip) at the edges, we can iterate the result Eq. ( 2.40| ) 
to any order: 

1(e) = (ie) n [ " dM (n) (ip)e^ . (2.41) 



The result Eq. (2.41) means that, when e — > 0, 1(e) = 0(e n ) for any integer n, i.e. 1(e) vanishes faster than any 
power of e. It does not mean that 1(e) is zero for any finite (but small) value of e. For instance, under stronger 
assumptions about the existence and properties of an analytic continuation of the function A(ip) in the complex tp 
plane it follows that /(e) ~ ie"? for some constants A and B. For reasonably small values of e such exponentially 
small contributions are numerically negligible. [We shall see later that the 'small parameter' e (or better e/B) is 
typically of order 1/(2ttN) where N = F 2 /F, is the instantaneous number of cycles.] 

We conclude, therefore that the integral / will be (in most relevant cases) numerically non-negligible only if the 
assumptions above are violated. In other words, the value of 1(e) will be dominated by the contributions coming 
from either (i) from stationary-phase points, ip(t s ) — 0, or (ii) from the edge points t a and/or tt>. Let us (for 
simplicity) assume that there is (at most) one stationary-phase point t s , and that it is of the normal parabolic type, 
i.e. tp(t) = ip s + ^ s (t - t s ) 2 + 0((t - t s ) 3 ) with ip s ^ 0. Let us also assume that a(t a ) ^ 0, a(t\,) ^ 0. [We 
maintain here for the moment, the assumption (iii) above about the regularity of the functions a(t), ip(t) in the 
closed interval [t a , t^]] Then, assuming analyticity of the involved functions, the mathematically most rigorous way of 
decomposing / as the sum (modulo nonperturbative small contributions of the type discussed above) of a stationary- 
point contribution /stationary and of edge contributions / e dge = /edge ^~ ^cdgc i s to deform the original (real) contour 
of integration into the complex plane [^5|,|lJ. The deformed contour must be such that near t s it leads to a basic 
integral of the type dx e~ bx [co + c\x + C2X 2 + ■••], while near each end point it leads to integrals of the type 
JJ dye~ cy [do + d\y + d2y 2 + ■••]. It is then easy to find the structure of the expansion of both /stationary and / c dgo 
in powers of e, as e — ► 0. For instance, it is convenient near t s to introduce a scaled variable: t — t s = e?r (before 
rotating t to complex values t = e ±1 ^x), so that the phase scales as 
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i> i>s , 1 V 2 , V s " , (3) 3 



(2.42) 



where ipg = d 3 tjjs/dt 3 . Expanding then the integrand of this r-integral in powers of s 2 leads to an integral of the 
symbolic type 



-^stationary ^ £ 2 J dx 6 



l + £ ia; odd + ex even + e*a; odd 



(2.43) 



where x odd (x cvcn ) denotes a sum of terms ~ x 2k+1 (x 2k ) . This yields (using the fact that the terms f_ a /^ 
are exponentially small) an expansion of the type 



dx < 



odd 



'stationary 



£ 5 [C + Cie + C 2 £ 2 + ■ ■ •] 



(2.44) 



The structure of the 'edge' contribution / c d gc can be obtained by a similar technique. Now the appropriate scaling is 



different, e.g. t — t a = er, and one ends up with integrals of the type dye v y n . This yields, 



/edge ~ £ [-Do + D l£ + D 2 e 2 ■ ■ •] 



(2.45) 



The aim of this preliminary discussion was to point out the s truct ures Eqs. (2.44) and ( 2.45| ) of the two main 
contributions to a generic oscillatory integral of the form Eq. ( |2.38 ). Note that while t he le ading contribution is 
given by the lowest order term in the stationary-point or sad dle-p oint expansion Eq. (2.44), the next-to-leading 
contribution comes from the lowest order edge-correction Eq. (2.45). One generally expects that / c d go will be only 
y/E, i.e. l/y2wN, sm aller than /statio nary We shall give below the explicit expressions for the first two ter ms in 

both expa nsions Eqs. (2.44) and ( 2.45| ). We shall see that each coefficient Co, C\, C 2 Dq, D\, ■ ■ ■ in Eqs 

and (2.45) is a combination of derivatives (of increasing total order) of ait) and ip(t) evaluated at t s for Eq. 
and at t a or tf, for Eq. ( 2.45| ). We n ote i n adv ance that, for actual calculations, the simplest way to evaluate the 



2.44) 



2.44) 



explicit forms of the expans ions E qs. ( 2.44 ) and ( |2.45| ) is not necessarily to follow the complex-contour route sketched 
above. In the case of Eq. (2.44) one can deal directly with the original stationary-point-expanded integral written 

as oc £5 f dr e l ^ sT [to + c\t + • • •], and in th e case of E q. ( 2.45| ) the simplest is to keep the boundary terms in 
the integration-by-parts approach Eqs. (2.40) and (2.41) given above for the simple case where they were neither 



stationary phase points, nor boundary contributions. 

To put in context the anal ysis t hat w e sha ll perform below, let us finally mention two serious limitations of the 
assumptions leading to Eqs. (2.44) and (2.45). First, in the analysis above, based on the introduction of the formal 
parameter e — > 0, we have assumed that the stationary- phase point t s was parametrically separated from the edges t a 

or tb, i.e. that \t s — t a \ and \t s — tb\ were much larger than the characteristic Gaussian width At = J E/ip s = O(^fe) 
associated to the stationary point. As we shall see, this limitation is unacceptable for the application we have in mind 
and we shall have to introduce new tools to overcome it. A second limitation (which compounds with the first and 
will lead to an unavoidable complexity of our treatment) is the seemingly innocent assumption (iii) above, namely the 
hypothesis that the functions a(t) and ip(t) are infinitely differentiable within the closed interval [t a , tb] (i.e. including 
also the end points). As we shall see, in the physically relevant case of a relativistic (adiabatic) chirp the functions 
a(t) and ip(t) will not be C°° at the physically imposed upper cutoff tb = £lso- This will require us to introduce new 
types of expansions and new tools to deal with the relativistic edge contribution in addition to the modification of 
the stationary-phase approximation needed in the case where t s is near the edge, in the sense that \t s — tb\ = O(yfe). 



III. IMPROVED STATIONARY PHASE APPROXIMATION FOR TIME- WINDOWED 

NEWTONIAN-LIKE SIGNALS 

A. The usual stationary phase approximation 



Let us begin by a quick recall of the usual treatment of the stationary phase approximation to a chirp. Consider a 
signal, 



where, 



h(t) = 2a(t) cos <f>{t) = a(t)e- l0(t) + a(t)e i ^ t '> 
d4>(t) _ 



dt 



2TrF(t) > . 



(3.1a) 
(3.1b) 
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We shall say that a signal is "Newtonian-like" if the instantaneous frequency F(t), defined by Eq. (3.1b), increases 
without limit when t runs over its full (mathematically allowed) range. (Note that we conventionally consider only 
positive instantaneous frequencies.) For instance, at the Newtonian order, the explicit forms for the chirp amplitude, 
phase and frequency of the gravitational waves are respectively given by: 



a(t) =C M {irMF(t)) 2 / 3 



4(t) = <Pc 



wMF(t) = 



(t e ~ t) 

5M 



5/8 



5M 



256(i c - t) 



3/8 



(3.2a) 
(3.2b) 

(3.2c) 



where Ai is the chirp mass given by M = rf/^m in terms of the total mass m and the symmetric mass ratio r\\ 
<j) c the gravitational wave phase at instant of coalescence t c and Cm is the product of a function of different angles, 
characterising the relative orientations of the binary and the detector, with the ratio A4/d where d is the distance to 
the source (see below). Note that the function F(t) increases without limit as t tends to the formal coalescence time 

Coming back to a general signal, the Fourier transform is defined by Eq. (1.1). Because the signal hit) is real, 
we have the identity h(—f) = (h(f))*. It therefore suffices to compute the Fourier transform for positive values of 
the frequency /. [Note that we use a lower case letter to distinguish the F ourie r variable / from the instantaneous 



frequency F(t).] The Fourier transform of a generic signal of the form Eq. (3.1a) reads 



h(f) = h.(f)+h+(f), 
where, = / dt a(t)e l(27r/t -* (t)) , 



M/) 



(3.3a) 
(3.3b) 

(3.3c) 



The integrands of h±(f) are violently oscillating and thus their dominant contributions come from the vicinity of the 
stationary points of their phase (when such points exist). When / > (which we shall henceforth assume), only the 
h—(f) term has such a saddle-point. Therefore, we can write the approximation 



ft(/)~/L(/)~ / dt a{t) e***® 

J — oo 

where , ipf(t) = 2ir ft - <j>{t) . 



(3.4a) 
(3.4b) 



The saddle-point of the phase ipf(t) is the value, say tf, of the time variable t where dxjjj (t)/dt = 0, i.e. it is the 
solution of the equation F(tt) = /. The dominant contribution to the integral Eq. (3.4a) now comes from a time 



interval near t = tf. When the second time derivative of the phase at the saddle-point does not vanish, i.e. when 



F(tf) 0, one can estimate Eq. (3.4a) by replacing ipf(t) [and a(t)] by truncated Taylor expansions near t = tf, 
namely 



a(t) ~ a(t f ) . 



(3.5a) 
(3.5b) 



(The zeroth-order term in Eq. (3.5b) is enough because the first order term a{tf){t — tf) vanishes after integration 
over t). This leads to a Gaussian integral 



Kf) 



dt a(t f ) e ^/(*/)-*^(*/)(*-*/) a 



(3.6) 



Evaluating this Gaussian integral, one finally obtains the well-known expression for the usual SPA (hereafter abbre- 
viated as uSPA): 



h uspa {f) 



a(t f ) 



3 »hM*/)-T/4] 



nt f ) 



(3.7) 
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where tpf(tf) is the value of ipf(t) at t = tf. 

The conditions for the validity of the SPA are usually assumed to be e% <C 1, £2 "C 1, where 



£1 



a{t)4>(t) 



£2 



1 Ht) 
2tt F 2 (t) 



2ttN 



(3.8) 



One can assess in a mor e precise quantitative manner the accuracy of the SPA by computing the leading correction 
t o the integ ral E q. (3.6). This leading correction will be given by keeping more terms in the Taylor expansions Eqs. 
(3.5a) and ( |3.5bD. To keep track of what one means by the 'next order term' in the SPA expansion it is convenient 
(as in Section ]!! D| above) to formalize the fast variation of the phases <f>(t) and ip(t) by considering an integral of the 
form I = J dta(t) exp(itp(i) / e) with a 'small' parameter e (set to unity at the end of the calculation). It is then easy 
to see [e.g. after the introduction of a rescaled time- variable: t — tf = e 1 / 2 r where tf denotes, as above, the saddle 



point of the p hase ipf{t)] that the leading correction to the result Eq. (3.7) will be of fractional ord er e [as ex hibite d 
in Eq. (2.44)] and will come from keeping two more terms in both the Taylor expansions Eqs. (3.5a) and (3.5t). 
Expanding in powers of e 1 / 2 leads to integrals of the form drr n exp(—inF(tf) r 2 ) with n < 6. Finally one finds 
that the sum of the usual SPA and of its leading correction is equivalent to multiplying Eq. (3.7) by the correcting 
phase factor e tS where (F^ denoting d 3 F/dt 3 ) 



8 = 



2nF(t f ) 



I'd 
2a 



la F 
2a F 




8 F 



(3.9) 



Therefore, a quantitatively precise criterion for the local validity of the SPA is £i oc = \S\ <C 1. In the case of power-law 
chirps, £i oc is equal to one- fifth of the criterion explicitly given in the recent study []30f| of the validity of the SPA. In 
the case of Newtonian chirps, Eq. (3.9) yields 



23/1 F \ _ 23 / 1 
~ 24 I ^F 2 j ~ 24 \9nN 



(3.10) 



Written in terms of v = (nmF) 1 / 3 this reads 8 = (92/45)t7w 5 which agrees with the corresponding result in [j3~L| . 
It is interesting to note that Eq. ( |3.10|) formally predicts, at the LSO, &lso = (4r?) x 0.58% which is quite small. 
Alternatively, one can say that Eq. ( |3.1C| ) predicts that even if the instantaneous number of cycles were as small as 
N ~ 1, the local correction to the SPA would be small (8 — 0.0339/7V). This result does not mean, however, that we 
can use the usual SPA Eq. (|]?]) to estimate with sufficient accuracy the FT of a real inspiral signal. Indeed, even if 
we w ere considering a Newtonian- like signal (where N stays away from zero at the LSO) the correcting phase factor 
Eq. (3.9) represents just the local correction to the SPA, i.e. the correction due to higher-terms in the local expansion 
near the saddle-point. There are also global corrections to the SPA coming from the entire integration domain, and, 
most importantly (as emphasized in |30j) from the end-points of the time-integration. In addition, there is also a 
correction com ing from the neglected contribution h+(f) in Eq. (|3.3b| ). Before considering them in detail, let us also 
note that Eq. ( |3.9| ) indicates that 8 blows up to infinity, at the LSO, in the case of a relativistic GW chirp (because 
F(t) ~ -Flso — a(^LSO — t) 1 !" 1 there; see below). This shows again that, independently of the problems linked to 



the time-windowing, relativistic signals will pose special difficulties, 
time-truncated Newtonian-like signals, by which we mean that N e 
cut-off. 



But let us start by studying the simpler case of 
= F /F stays away from zero at the upper time 



B. Beyond the usual stationary phase approximation 

We therefore consider time-domain signals of the form 

h(t) = 2a{t) cos 4>(t)B(t ma x - t) , (3.11) 

where 8 denotes the Heaviside step function, (9(x) — 1, for x > and 6(x) = 0, for x < 0). This time-windowing has 
three effects: (i) it induces oscillations in the usually considered frequency-domain / < F max , (ii) it generates a tail in 
the usually disregarded frequency-domain / > f max , and (iii) it renders non- negligible the 'non-resonant' contribution 
h+(f). Here, F m ax denotes the instantaneous gravitational wave frequency reached at t = t 

max, i.e. fmax — F'(^max). 

The main purpose of the present paper is to model and estimate, analytically, as accurately as possible all these 
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effects. The case where the saddle-point tf is (below and) far away from the upper-cutoff i max has been recently 
considered in jjlj. However, this case is not the physically relevant one. As pointed out in DIS, the case where the 
usual SPA becomes unacceptably inaccurate is the case of massive systems for which the signal emits not many cycles 
in the detector's bandwidth before crossing the last stable orbit. In this case the most important frequencies are 
located around the effective-cutoff frequency F max ~ -Flso (and as we shall see below, it is important to estimate the 
Fourier transform accurately, both for / < F majt and for / > F max ). Here, we shall provide an approximate analytical 
treatment valid in this crucial range of frequencies. 

Let us first state clearly our notation. We decompose the FT of the time- windowed signal Eq. ( |3.11 ) as 



M/) = M/) + M/), 



where, h-(f) 



h+(f) 



«/> f (*) 



— oo 



dt a(t) e 



dta{t) e^/ (t) 



TPJ{t) = iP f (t) = 2Trft-<f>(t), 



(3.12a) 
(3.12b) 

(3.12c) 

(3.12d) 
(3.12e) 



We shall refer to the contribution h-(f) as the 'resonant' contribution because the equation defining the saddle point it 
contains, F(tf) — /, corresponds to a resonance between the Fourier-frequency / and the instantaneous gravitational 
wave frequency F(t). 



C. Edge contribution to the non-resonant integral h+(f) 

Before dealing with the oscillatory and tail corrections to the resonant contribution h-(f) of h(f), we shall first 
deal with the non-resonant contribution h + (f). When / > 0, the phase i/j~f(t) in contrast to ipj(t) has no stationary 
point. Let us therefore consider the general problem of approximating an integral of the form 



I 



dta{t)e^ {t) 



(3.13) 



where ip(t) is always monotonically varying (ip ^ 0). We can then use tp as integration variable. This yields, as above, 



1 = 



/ dip A(ipy^ , 



(3.14) 



where A(ip) = a(t)/ip(t) where t(ip) denotes the (unique) solution in t of if) = ip(t) and where tb a = ib(t a ), t/Jb = 

Jt(V-) 

ip(tb). One can treat ip as a fast varying phase, so that by comparison, the amplitude A(ij)) varies slowly when ip 
varies by 2-7T (as above this could be formalized by the formal replacement ip — > ip/s). In other words, instead of a 
SPA, we are in a WKB like approximation where one can expand in the slowness of variation of A(tp), i.e. we can 
expand in successive derivatives d n A(ip) / dip n . This exp ansion is obtained by successively integrating Eq. (2.39) by 
parts [using e 1 ^ = /i)]. Contrary to Section II D above we now keep the edge contribution coming from the 

boundaries. For instance at second order this leads to 



I = 



■—[AW + iA'M] 



+ i 2 dil>A"($)e i + . 

a •'i'a 



(3.15) 



It is easy to re-express this result in terms of the original time variable t by replacing A(ip) = a(t)/ip(t) and d/dip = 
(ijj)~ 1 d/dt. We deduce from Eq. ( |3 .15 ) the full 'edge' contribution to the integral I (as explained in Section |ll D| , the 
'bulk' contribution is exponentially small): 



-'edge — 



(A(V) + iA'(if>) + ■■■ + i n A^ (V) + ■ • •) 



(3.16) 
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This is the explicit form of the parametric expansion sketched in Eq. (2.45). It can be expressed in terms of the 
time derivatives of a(t) and ip(t) by using the replacement rules just mentioned. In particular, at the leading and 
next-to-leading order it reads explicitly 



ii>(t) 



1 + 



1 U{t) a{t) 
i#) \i>{t) o(t) 



If 



we apply this general result to h + (f) Eq. ( 3.12c ), we get the estimate 



h + (f) ~ ~ "f max) e ^(W) 

^(^max) 



i + 



2™(/ + F max ) 



1 / g(W) _ dft max ) ' 

#j(imax) \^/(<max) a (W) j 
1 / ^inax ^(^max) \ 

27ri(/ + F max ) I (/ + F max ) a(t 

max/ / 



(3.17) 



(3.18a) 
(3.18b) 



Note that, in the approximation where the amplitude is taken as Newtonian-like, i.e., a(t) — Cvj? oc (F(t)) 2 ^ 3 , the 
last term in Eq. fl3.18b|) reads a(t max ) / a(t max ) = § F max /F max . 



D. Improved stationary phase approximation when / < 



Let us now consider the dominant contribution ft,_(/), Eq. (3.12b). We st art by considering the case where the 
Fourier variable / [of the Fourier transform of the time- windowed signal Eq. ( |3.11 )1 is smaller than -F max , but near 
F max . As will become clear from the formulas below, the interval around _F max where it is needed to improve the usual 
SPA, is the range, 



1/ - F n 



< 



(few) y/I^tf). 



(3.19) 



In the case where / is in the range Eq. ( |3.19| ) with / < F max , there is a sad dle-p oint tf in the first term of the exact 
integral Eq. ( 3.31 ), and o ne can still use the parabolic approximation Eq. ( [3.5a ) to the phase ipf(t), and the lowest 
approximation Eq. (3.5b) to th e am plitude a(t). [Indeed, the work of the previous Sect ion h as sho wn t hat the local 
corrections to the integral Eq. (3.4a), coming from the inclusion of more terms in Eqs. ( |3.5a| ) and ( |3.5b| ), were quite 
small as long as N > 1]. Therefore, in this case, the resonant contribution to the Fourier transform becomes: 



h-(f) 



alt a(t) e"^ (t) , 



a(t f )e 



dt e 



-inF{tf){t-t,f 



(3.20a) 
(3.20b) 



The crucial difference between the Eqs. ( |3.6| ) and (3.20b) is that the full Gaussian integral has become a com- 
plex Fresnel integral which may be evaluated in terms of the complementary error function. Let us recall that the 
complementary error function erfc(z) = 1 — erf(z) is defined by 



2 r+oo 

erfc(z) = —= I i 

V 7T 



dx . 



(3.21) 



It takes on the real axis the particular values erfc(+oo) = 0, erfc(0) = 1, and erfc(— oo) = 2. By rotating the 
integration contour in the complex plane (x = — e~£) and shifting the new integration variable £, we get the following 
useful integration formula: 



dj . e -i(ar 2 +2br+c) 



1 n 



e 4 e 



erfc 



-e" y/a(T m + -) 



The formula Eq. (3.22) motivates us to define the following auxiliary function 

C(C)^erfc(e^c) ■ 



(3.22) 



(3.23) 
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It is useful to note that C(+oo) = 0, C(0) 
expansions of C(C) as (, — > ±00 are 



1/2, and C(— 00) = 1. Moreover, the leading terms in the asymptotic 



C->-oo, C(C)~1 + 



C 



-co, 



C(C) 



2y^ C 

*C 2 



1 e 



(3.24a) 



(3.24b) 



2^ C 

As the auxiliary function C(£) plays an important role in our work we plot it in Fig. g. From the Figure we note 
that the real part of the above complex function is a softened version of a step function 9(—() and the imaginary part 
an oscillating function vanishing in the limits ±00, as well as at £ = 0. 



Armed with this definition one finds that the right hand side of Eq. (3.20b) yields the approximation 

M/) ^c(Co(/))^ uspa (/), 



where, 



Co(/)= \I^F{tf){t f -t max ) 



(3.25a) 
(3.25b) 



In words, when / < F max on e ca n correct for the 'edge effects' caused by the cutoff at i max by mul tiplyin g the usual 
SPA A uspa (/) given in Eq. (3/7) by a complex 'correction factor' C(Co(/))- The expression Eq. ( 3.25a ) gives very 
good overlaps with the exact discrete Fourier transform (DFT) of the time-wind owed s ignal Eq. (3.11). However, it 
is possible to do even better by a slight modification of the argument Co(/)j Eq. ( |3.25b| ). 

To understand how a slight modification of the argument £(/) of the auxiliary function C(C) can improve both the 
visua l agreement (even quite far away from F ma x) and the overlap wi th the exact DFT of the time- windowed signal Eq. 
( 3.11 ) we have to take into account the asymptotic expansion Eq. (|3.24a ). Indeed, on the one hand, when inserting 
the expansion Eq. ( 3.24a| ) into Eq. ( 3.25a ), using Eq. ( |3.7j ) for h uspa, (f) and allowing for a more general frequency 
dependent argument C(/)> we nn d that h(f) differs from /i uspa (/) (in the domain £(/) — > —00 i.e. f <§; F max ) by 
a correction term proportional to e~™/ 2 e l ^ f ^ tf )~^ 2 1 /(*. On the other hand, a different way of estimating this edge 
correction consists in writing Eq. ( |3.20b ) as an integral between —00 and +00 minus a "correcting" integral between 
^max and +00. As was discussed in Section [II D| the latter integral can be estimated by successive integration by 
parts. This gives [see Eq. (3.17)] a first order correction term proportional to e + I71 7 2 q(t max )e^ / ^ max ^ /ipf (£max )- The 
phasing of this edge correction can be made to agree perfectly with the phasing predicted by the form Eq. (3.25a) 
written with a generalized argument ((/) if (ipf(tf) — ( 2 ) = i/jf(t max ). This leads us to define, in the domain C < 0, 
i.e. f < -Fmax, the new argument 



-\Ji>f(tf)-i>f(tr. 



(3.26) 



In the left part o f the crucial region E> 
Co(/) Eq. (3.25b), as is seen from Eq 



{3T9) the argument C<(/) Eq. (3.26) is nearly identical to the previous result 
5a). However, we have checked that the replacement of Co by £< improves both 



the visual agreement and the overlap with the exact h(f). Let us also note in passing that an amplitude proportional 
to of the correction term to h uspa (f) derived from the asymptotic expansion of Eq. (3.24a) is consistent with 



the different analytical treatment used in |31j which was valid only for / <C F max , i.e. large, negative £. By contrast 
our approach based on the function C(£) is adequate in the full range — 00 < ( < without exhibiting any fictitious 
blowup at C = (remember C(0) = 1/2). [Our approach is also valid in the region £ > 0, i.e. f > F max , but for an 
improved treatment of this domain, we shall find it convenient to modify further the argument C(/) m the following 
Section.] 

Summarizing: we propose as final result for the resonant part of our improved SPA for Newtonian-like signals (or 
inSPA) 



/ < F n - 



h™r=c«<(f)) 



F(t f ) 



(3.27) 



The corresponding total improved approximation ft, lntot to the Fourier transform h = h-(f) + h + (f) is the sum of 

Eqs. ( 3. 18b| ) and ( 3.27 ). Note that the ratio fa _/fa + is (w hen t f is near £ max ) of order 47rF max / \J F ma x = 47r v / iV r ma x 
(this is consistent with the e scaling of Eqs. (2.44) and ( 2.45| ), remembering that e ~ 1/2-7^). The contribution 



of h + is expected to be non-negligible only for signals which are really discontinuous in time. As the real signal 
(whatever be the subsequent plunge signal) will be continuous (and even smooth) it is clear that one should not add 
any contribution from h+ when applying our above treatment to real signals. In fact, we shall see below that, even 
for discontinuous signals, the addition of h + has only a minute effect on overlaps. 
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E. Approximate Fourier transform when / > i 7 , 



Let us now consider the evaluation of h-(f) in the case when the Fourier variable / is larger than F m:ix [but near 
F max , in the sense of Eq. ( 3.19] )]. In that case the integral Eq. (3.12b) giving h-(f) no longer has a saddle-point. 



However, it 'nearly' has a saddle-point and therefore we expect that 



h-(f) 



dta(t) e^f {t) 



(3.28) 



will still dominate over h+(f). One could think of two ways of analytically approximating the integral Eq. ( |3.28 ). A 
first way is to still use the fact that (for Newtonian-like signals where the mathematical function F(t) continues to exist 
and increase beyond t = t max ) though there is no saddle-point in the domain of integration [— oo,t max ], there exists 
a nearby saddle-point of the analytically continued phase function ijjf(t). More precisely, for Newtonian- like signals 
the mathematical equation F(tt) = / still defines a unique value tt (with tt > i max when / > -F max ). Capi talizing 
on the existence of this nearby saddle-point one can still try to insert the expansions Eqs. (3.5a) and (3.5b). This 
leads to a result of the for m Eq. (|3.27| ) with the correction factor Eq. ( 3.25b| ) , i. e. now considered for positive values 

of the argument Co(/) = \J ^F(tf)(tf ~ tmax)- In other words, a simple uniform expansion to h—(f) on both sides of 
/ ~ Fmax would seem to be simply 



\fHtf) 

with Co(/) = \JnF(t f )(tf - i max ) . 



(3.29a) 



(3.29b) 



Here 'cspa' means (zeroth order) corrected SPA. Note that when / > F max , tt, and therefore all the quantities 
evaluated at tf, are defined by using the (supposedly existing) analytic continuation of the mathematical function 
F(t) beyond t = i max . 

In our firs t attem pts at improving the SPA in presense of a time-windowing we came up with the simple proposal 
Eqs. (|3.29aD - (|3.29b|) and it gave excellent overlaps with the exact Fourier transform. However, we realised later that 
we could further improve on th is simple p ropo sal. We already stated that for / < F max our best proposal is to 
modify the argument Eq. ( |3.29b|) into Eq. ( 3.26 ). In the case where / > -F max our best proposal is neither to use the 
straightforward argument Eq. ( |3.2'9"b" ), nor the "improved-phasing" ar gumen t Eq. ( p.26[ ) with a positive sign in front 
of the square-root [which, however, still improves over the choice Eq. ( |3.29b )] but to follow a different tack which will 
turn out to be useful when considering the case of relativistic-like signals in the next Section. 

To motivate our proposal in the case / > -Fmax, let us remark that the integral to be ap proxi mated, i.e. Eq. 
( |3.28 ) having no saddle-point in the domain of integration, is formally of the general type Eq. ( |3. 14 ) with the phase 
t Mt} — jifit ) bein g a monotonically increasing function of t. The important information we wish to deduce from Eqs. 
(BT5|) and ( jug ) is that there exists an expansion (valid when / ^> F max , i.e. Co(/) is large and positive) in which 
h—(f) is entirely expressed in terms of the values of the functions ipf(t) and a(t), and their derivatives, evaluated at 
the edge point t = i max . This contrasts with the 'corrected' result Eqs. ( 3.29a| )-( |3.29b| ) which relied on the existence 
of the functions ipf{t) and a(t) in the "unphysical" region t > t max . This motivates us to look for an approximation 
to Eq. (2.39) val id al l over the d omain Co(/) > [and not only when (o ^> 1 which will be seen to be the domain of 
validity of Eqs. ( |3.15| ) and (3.17)] but expressed entirely in terms of the edge valu es of tpf{t) and a(t). We propose to 
define such an approximation by replacing the phase and amplitude in Eq. ( 3.28 ) by 



1pf(t) ~ ^>/(tmax) + 27r(/ - F max )(t - t max ) - TtF(t max )(t - t max ) 2 

a(t) ~ <z(£ max ) . 



(3.30a) 
(3.30b) 



Thanks to the parabolic nature of the approximation Eq. (3.30a) this again leads to an incomplete complex Gaussian 
integral (i.e. a Fresnel integral) which can be evaluated as before in terms of the complementary error function. Using 



Eq. (3.22), this leads to our final proposal for the (nearly) resonant part of h-(f) for Newtonian-like signals 



f>F n 



^(^max) 



c pa (/)^(c>(/)) 

sJF{t 

^ max) 







exp i 



Af - -F max ) s 

^(^max) 



tt/4 



C>(/) 



^(^max) 



(3.31a) 
(3.31b) 
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Note that, in the parabolic approximat ion wh ere Eq. ( |3.5a ) or Eq. (3.30a) hold, the f uncti on C>(/) i s a P~ 
proximately equal both to Co(/)> Eq. ( 3.29b ), and to the analytic continuation of Eq. ( 3.26 ), i.e. C<(f) = 
sign(/ — F max ) y/ 'tp f (t /) — ^/(W)- Note also that the phase factor in Eq. ( |3. 31a ) (which is explicitly expressed 
in te rms o f edge quantities) is nearly equal to the analytic continuation of the usua l SPA phase factor appearing in 
Eq. (3.27), i.e. exp [iipf(tf) — wr/4]. Finally, as required, the expressions, Eqs. (3.27) and ( 3.31a| ) match continuously 
at / = -Fmax with common value 



(-^inax ) 



, mspa/ 

< \^ max } 



0* (Wx) f (t max ) -7T /4) 



^(^max) 



(3.32) 



Summarizing: our best analytical estimate for the Fourier transform of discontinuous Newtonian-like signals is the 
sum 



t(/)=/i mspa (/) + /i od gC(/)5 



(3.33) 



where /i!j_ se (/) is approximated by Eq. (3.18b) and where ft,'" spa (/) is given, when / < F max by Eq. (3.27) and for 



/ > -Fmax by Eq. ( 3.31a ). As stated earlier, we shall in fact recommend that the edge correction h e + &e be not included 



when applying our result to real signals (we shall also see that it brings only a negligible improvement to the overlaps 
of time- windowed signals). 



F. Comparison between the improved SPA, the usual SPA and the 'exact' SPA (numerical DFT) 



Before proceeding to a quantitative comparison of the various approximants in Table [n|, it is important to remark 
that one needs to be more specific when using the terminology uSPA. One could compute the uSPA truncated at the 
-Flso that we refer to as uSPAw ( where 'w' stands for 'windowed') or the uSPA truncated at the Nyquist frequency 
that we designate as uSPAn (where 'n' stands for 'Nyquist'). In Table || we have listed the overlaps, as defined by Eq. 
(2.8), of a signal model generated in the time-domain and then Fourier transformed using a numerical DFT algorithm^ 
with the same signal model but directly generated in the frequency-domain using the usual SPA (uSPA) , the corrected 
SPA (cSPA) and the improved Newtonian SPA (inSPA) discussed earlier. For simplicity, we consider only equal-mass 
systems (r) = 1/4) and parametrize them by the total mass m = mi + m,2- The total mass is the crucial parameter 
which measures the location of the Flso with respect to the bandwidth of the detector. The parameter r\ is also 
important because it determines the number of cycles near the LSO [Eq. (2.27) shows that N(F^so) scales as l/rj\. 
The worst case (for the sensibility to the shutting off of the signal after the LSO) is r\ = >7max = 1/4, and this is 
why we focus on this case. [We are also motivated by the fact that 1/4 being the maximum value of the function 
<7(mi,m2), the observed values of ?7, corresponding to a random sample of of mi and m2, are expected to have an 
accumulation point at ?7max = 1/4.] As we have checked, if our filters exhibit good overlaps for r\ — 1/4 they will have 
even better overlaps for 77 < 1/4 and the same value for m. 

The error function needed in computing C(f) is numerically computed using the NAG library S15DDF. The overlaps 
are shown for the usual SPA with a frequency-windowing (uSPAw) together with the overlaps for the uSPAn, cSPA 
and inSPA, computed up to -Fkyquist- Table || shows that the improvements on the SPA that we propose in this paper 
(both the simple cspa and our final inspa) succeed very well in modelling the edge effects due to time- windowing. The 
overlaps in the case of cSPA/inSPA are better than 0.99 for equal-mass systems with total mass m < 40M Q . For a 
system of m = 40Af Q uSPAw gives an overlap of 0.8589 resulting in a loss in the number of events by 37%. Although 
the overlaps of cSPA and inSPA seem to be always about the same, we think that inSPA is a better representation of 
h(f)] it has better overlaps in the case of the most massive systems (see the first lines of Table ||), and, as shown by 
Fig. ||, it captures better the decay of h(f) beyond F ma x- In the Table for the usual SPA we have listed the overlaps for 
uSPAw i.e., uSPA terminated at / = Flso = 4400(m©)~ 1 Hz and the overlap (uSPAn) up to the Nyquist frequency 
/Nyquist = 2 kHz. As is very clear from these entries, windowing of the SPA improves the overlaps for massive systems 
very much. As remarked earlier, this is why in DIS, while comparing the DFT to the SPA, the uSPAw was used. On 
the other hand, computing overlaps up to the Nyquist frequency, i.e. uSPAn, produces much smaller overlaps. 



6 This defines for us the 'exact' Fourier representation of the signal, after due care has been taken to use a smooth time- window 
below f s , and a high enough sampling rate. 
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To understand this further we plot in Fig. || the power per logarithmic bin of the squared SNR, dp 2 /dlogf — 
/|/i(/)| 2 /S , T1 (/), which is the Fourier-domain quantity of most significance when discussing overlaps. We compare this 
quantity for various approximations to the Fourier-transform of an (arbitrarily-normalized) time-windowed signal: 
DFT, uSPA, cSPA and inSPA. In the important range of frequencies our best analytical approximant inSPA agrees 
with the exact result (FFT) quite well, the uSPA grossly overestimates and cSPA somewhat underestimates the 
actual signal power. This is why, though analytically continued up to /Nyquist, the uSPAn returns a smaller overlap 
as compared to the windowed SPA (uSPAw) because it overestimates the power in the signal beyond FLso [M • 

In all the comparisons above, it is worth stressing that the FFT calculation is delicate: The 'exact' time-domain 
chirp contains an infinite number of cycles in the far past, with instantaneous frequencies tending to zero. As what 
happens to frequencies below the seismic cut-off / s = 40 Hz is not physically important, we wish to simplify the 
numerical calculation of the FFT by essentially discarding the (infinite) part of the signal, having instantaneous 
frequencies F(t) < f s = F(t s ). We started doing that by simply time-windowing the signal for t < i min < t s by a 
sharp, lower time-window 6(t — i m in)- However, this method introduces physically spurious oscillations (which are the 
lower-cutoff analogue of the physically important upper-cutoff oscillations) present in both h+(f) and h-(f). One 
way to deal with this problem is to subtract out from the FFT these spurious edge oscillations by using the general 
formula Eq. ( |3.17 ), which in the present context, can be applied both to h^ FT (f) and hZ FT (f). For instance, to 



lowest order the FFT corrected for these oscillations would read 



ejected = ~h FFT + A+ in + A min , (3.34a) 
where, A± n = — ^ j+t^ , (3.34b) 



where Finm = ^(imin) < fs [a better approximation can be derived from Eq. (3.17)]. However, it seems better to use 
an alternative approach which does not require one to correct by hand the FFT. The alternative approach we have 
actually used in our calculations consists in imposing a smooth (rather than a sharp) lower time- window on the exact 
chirp, acting below t s and in a smooth-enough manner that it does not introduce spurious edge-oscillations in the 
frequency-domain. The smooth time-window that we used consists in multiplying the chirp by the function 

e z + 1 t — 1\ t — t% 

which smoothly interpolates between when t = t\ + and 1 when t — t2 — 0. We used t\ such that F\ = F(ti) = 30 
Hz and ti = t s i.e. Ffa) — f s - Moreover, we need to be careful with sampling and phase factors to correctly reproduce 
the edge correction to ft.^ ge . 

In addition to the comparative evaluation of the various approximants, Table || also provides a numerical proof 
regarding the effect of h c + 6C on the overlaps. It is quite important to note that the inclusion of the non-resonant edge 

term h 6 ^^ has only a very minute (but positive) effect on overlaps. This is good news for our formal time-windowing 
ansatz, because we expect that this contribution will be (exponentially) negligible in the case of real (continuous) 
signals. We interpret the fact that even for our formal discontinuous model h e + ge is negligible^ as a confirmation 
that our improved SPA can adequately model not only signals that vanish after the LSO, but also signals that shut 
off rather quickly (on the F£g time-scale) after the LSO. It leads us also to propose, finally, to use as analytical 
representative of the FT of real signals the ft," lspa part of our formula above (without the edge term). [To simplify the 
notation, we shall henceforth drop the extra subscript minus on /i mspa .] 

In computing the above overlaps we have matched all the parameters of the two waveforms, including the time 
of arrival and the starting phase^. The overlaps in this Table as well as all other Tables in this paper are found to 
be insensitive to the sampling rate at the level of a fraction of a percent, provided that it is large enough to obey 
Shannon's sampling theorem. 



7 Note that our statement here is only that ft.? dse (/) can be effectively omitted without significantly worsening the overlaps. 
We are not claiming that /i+ ge (/) is pointwise numerically negligible compared to /i_(/). Indeed, because the instantaneous 
number of cycles is rather small near the LSO, our analytical estimates above show that /i! dge (/) is not very much smaller than 
h-(f) near / = Flso- 

8 The lag is set equal to zero in testing the accuracy of the Fourier representation but chosen optimally when testing faithfulness 
of a family of templates e.g. in Sec. M. 
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IV. IMPROVED STATIONARY PHASE APPROXIMATION FOR RELATIVISTIC SIGNALS IN THE 
ADIABATIC APPROXIMATION: THE SPP-APPROXIMANTS 



Though one might a priori think that it is a simple matter to generalize the improved SPA discussed above for 
Newtonian-like signals to the relativistic case, it does not turn out to be so. What complicates matters is that there 
are serious qualitative, "non-perturbative" differences between the two cases: first, the value of F(t) formally tends to 
+00 at the last stable orbit (LSO) which physically defines the upper-cutoff i max of the inspiral signal, and, second, 
the mathematical function F(t) does not admit a unique real analytic continuation beyond t m ax = ^lso- (These two 
facts are evidently related; indeed we shall see that F(t) behaves in the non-analytic manner F(t) ~ c\ +C2(£lso — t) 1 
when t — > £lso)- Remembering the crucial role of a finite F(t) in the results Eqs. (3.27) and (3.31a), it is clear that 
we need to tackle afresh the problem of finding a good, analytic approximation to h-(f). Similarly, in view of the 
appearance of F +1 (t max ) in the next-to-leading contribution to h+(f) Eq. (3.18b), we shall also need to revisit the 
calculation of h+(f) (though we shall, again, find that it makes only a negligible contribution to the overlaps.) 



A. The phasing formula for relativistic signals in the adiabatic approximation 

To extend the treatment of the previous Section and go beyond the Newtonian approximation, let us begin with the 
phasing formulas for gravitational waves from compact binaries written in a parametric form in terms of the variable 
vp = (ttitiF) 1 / 3 defined by the total mass m = mi + m 2 and instantaneous gravitational wave frequency F 

t(v F )=t hS o+m dv—j-^, (4.1) 



E'(v) 



0lso + 2 / dvv 3 — ff, (4.2) 



where E(v) is the dimensionless energy function related to the total relativistic energy or Bondi mass by -Etotai 
m(l + E), T(v) the flux function denoting the gravitational wave luminosity of the s ystem an d £lso is the time and 



^lso is the phase of the signal when v — ulso- The parametric representaion Eqs. (4.1 (4.2 of the phasing formula 
<f> = 4>(t) holds under the assumption of 'adiabatic inspiral', i.e., that gravitational radiation damping can be treated 
as an adiabatic perturbation of a circular motion. See fdljfl for a treatment of radiation damping going beyond this 
approximation. 

In the restricted post-Newtonian approximation, one uses a Newtonian approximation for the amplitude [ p5[ . 
However, in order to extract an inspiral signal that may be buried in noisy data by the method of matched filtering, 
we need to employ post- Newtonian accurate representations for the two functions E'{v) and F{v) that appear in the 
above phasing formulas. To any approximant Ea(v), Ta(v), correspond [by replacing E(v) — > Ea(v), T(v) — > Ta(v) 



in Eqs. (f4.1|) and (4.2)] some approximate parametric representation t — tA(v F ), 4> — <Pa(vf), and therefore a 



corresponding approximate time-domain template 

h A = h A {t;C,t LSO ,(f)LSO,m,r 1 ), (4.3) 

obtained by replacing Vp, in the following ^-parametric representation of the waveform 

h A (v F ) =Cvp cos (/)a(vf) , (4.4) 

by the function of time Dp — VA{t) obtained by inverting t = tA(vp). 

The standard approximants for E(v) and !F(v) are simply their successive Taylor approximants Er n and Tt^ 
respectively. The DIS strategy for constructing new approximants to E(v) and T(v) is two-pronged: Starting from 
the more basic energy-type and flux-type functions, e(v) and l(v) jl3| we construct Pade-type approximants, say 
e P„ 1 ^p„: of the "basic" functions e(v), l(v^. We then compute the required energy and flux functions entering the 
phasing formula. The successive approximants E[ep n ] and 3~[ep„,lp„] have better convergence properties than their 



9 For explicit formulas representing E(v) and T{v) see Eqs. (3. 8), (4. 2) and (4.3) of DIS. The associated e(v) and l(v) functions 
are given by Eqs. (3.7), (3.9) and Eqs. (4.4)-(4.9) in DIS . See also Eqs. (3.5), (3.11) and Eqs. (3.18)-(3.23) there. 



24 



Taylor counterparts Ex^ler^] and J-t„ [er„ , lr n ] ■ In DIS we were working directly with the time-domain signal h(t). 
As explained above this necessarily requires a numerical inversion of the parametric representation t = t(vp). By 
contrast, if one wants to compute the usual stationary phase appro xima tion of h{t) — Cv 2 F (t) cos <fi(vF(t)) there is 
no need to invert this parametric representation. Indeed, from Eq. (3.7), it is sufficient to know the instantan eou s 
amplitude and the phase at the time tf where / = F(tf). This time is simply given by the same expression Eq. (4.1) 
as above with the replacement of Dp = (nmF) 1 ^ 3 by vj = (ttto/) 1 / 3 , i.e. the stationary point tf is given by 



^LSO 



E'(v) 



dv . 



One then substitutes this value of tf in Eq. (4.2) to compute the phase ^/(i/) 
component: 

rVhso E'(v) 

il> f {tf) = 27r/i LSO - 0lso + 2 (vj- v 3 )—ffdv . 

In terms of these quantities one has 

fi uspa (f) = -C — Qityfttf)-^] 
2 



(4.5) 

2irftf — (j>(tf) of the Fourier 

(4.6) 

(4.7) 



The inclusion of relativistic effects in h uspa (f) is then simply accomplished by using relativistic accurate expressions 



for E'(v) and 3-{v) in the formulas giving tpf{tf) and F(tf). The coefficient C, in Eq. (4.4), determining the actual 
amplitude of the waveform reads: 



C(r,z,0,0>) = (477) (J£)C(i,6,h$), 
where d is the distance to the source, and where 



C7(i,M,V) = Va 2 + b 2 , 



with, A = - (l + cos 2 i) F + ; B = cosi F x 



with the beam-pattern factors 



F+ (0,4>,ip) = - (f + cos 2 6) cos 20 cos 2ip - cos 9 sin 2cj) sin 2ip, 
F x (0, 4>, ip) = - (I + cos 2 9) cos 2<[> sin 2^> + cos 9 sin 20 cos 2-0. 



(4.8) 

(4.9a) 
(4.9b) 

(4.10a) 

(4.10b) 
(4.10c) 



In these formulas the angle i denotes the inclination of the orbit with respect to the plane of the sky, and the angles 
9, 0, and ip parametrize both the propagation direction and the polarization of the gravitational wave with respect 
to the detector (see || for exact definitions; we added a bar over and ij) to distinguish them from the GW phase 
and Fourier phase ip respectively). Performing averages over the angles in the squared SNR leads to : 



(El) 



',<j>,ip 



and finally 



(C 2 



(F 2 X 



\<f>,1p 



4 

25' 



1 

5' 



(4.11) 



(4.12) 



We are finally in a position to write down the rms and ideal SNRs. For a binary at a distance d from the earth consisting 
of stars of individual masses mi and m2 (total mass m = mi + ?7i2 and symmetric mass ratio 77 = mim^/m 2 ) the 
rms and ideal SNRs, obtained by using the rms and ideal values of C, namely C — 2/5 and C — 1, respectively, when 



replacing E q. (fO |) in Eq. (2.20), or equivalently, when replacing a(/) 
(with Eq. ( p.26[ ) and a truncation at .Flso), are given by 



(1/2)C?j 2 (/) = 2 V md- 1 Cv 2 (f) in Eq. ( gjgg ) 



25 



5/6 f ^ N 1/2 



dlT 2 / 3 V 15 







1/2 



5 

Pidcal = TTPrms- (4-13) 



Note that the SNR depends only on the combination M. = mrf/ 5 - the chirp mass (see e.g. bM), and that the first 



Eq. (4.13) is equivalent to Eq. (La). 



Let us next delineate the qualitative differences between the relativistic and non-relativistic cases by considering 



the function appearing as denominator in the uSPA, Eq. (4.7) 



[ > 2vr dt 2 ~ irm 2 E'(vY [ ' 

At the LSO, the gravitational wave flux T(v) is finite (it blows up only later, when reaching the light ring jl3|]) while, 
by definition, E'{v) vanishes linearly, E'(v) oc v — «lso- As we shall see below this means that F(t) blows up as 



(*lso — t) I . A consequence of this blow up is that the last two terms in Eq. (3.9) blow up like (iLSO — t) 



-3/2 

confirming the need for a special treatment of the Fourier transform near the LSO. We are here speaking of the exact 
behaviour of the functions E{v) and F{v), as supposedly known from combining the test-mass limit results fl39| with 
the best available results on the physics underlying the existence of the LSO (l3J, and the emission of gravitational 
waves in comparable mass systems fl24j . In DIS, we have incorporated this information so that all the P-approximants 
Ep n = E(ep n ), J-p n = J 7 [ep n , lp n ] that we define share, with the "exact" functions E and T the crucial properties 
mentioned above (i.e. finite ^F(wlso) and E'{v) aw - ^lso)- The (less-convergent) successive T-approximants Et„ 
and Tt„ do not incorporate this information exactly, and only as n increases they tend to incorporate it. In our 
opinion the T n approximants disqualify as 'relativistic' approximants since they do not consistently incorporate the 
expectation (based on several different methods; see references injEcj) that the frequency at the LSO is (for any 



V < 1/4) numerically near the Schwarzschild-like prediction, Eq. ( 2.14 ). Indeed, if we define the 2PN Taylor estimate 
of -Flso by the value of v = (nrnF) 1 ^ where the straightforward Taylor approximant Et 4 (v) = X^t=o Ek(v) vk reaches 
a minimum, we find, e.g. that (i) when m — 40M Q and rj = 0, Fr 4 = 200 Hz, very different from the exact value of 110 
Hz, and (ii) when m — 40M Q and rj = 1/4, that F£g = 221.4 Hz, very different from the other predictions f£§ = 143 

Hz and F^q @ = 118.6 Hz. We compare and contrast in Fig. ^ the Newtonian and relativistic behaviours of the 
wave amplitude and instantaneous frequency F(t) during the last couple of orbits before the LSO. The blow up of 
F(t), i.e. the fact that the slope of F(t) becomes vertical is an effect which is localized in the last part of the last cycle 
before the LSO. Note also in Fig. [7] that a less localized consequence of this blow up is that the average frequency a 
few cycles before the LSO is smaller (for a given -FLso) m the relativistic case, than in the (unphysical) Newtonian 
one. Note that the physical origin of the blow up of F is that, just before the LSO the 'effective potential' for the 
radial motion becomes very flat (before having an inflection point at the LSO). In picturesque terms, the radial motion 
becomes "groundless" at the LSO. Evidently, the blow up of F is due to our use of the 'adiabatic' approximation 
down to the LSO. In reality, radiation reaction will cause a progressive transition between the inspiral and plunge 
which will modify the evolution of F(t) in the last cycle before the LSO. We shall discuss this issue in detail in a 
forthcoming paper 143] and subsequently its data analysis consequences. 



B. Edge contribution to the non-resonant relativistic h+(f) 



As in the Newtonian-like case we decompose h(f) in two contributions, Eqs. ( 3.121j ) and ( [3.12c ). The non-resonant 
contribution h+(f) will be dominated by the 'edge' contribution to an integ ral of our usual type Eq. (2.38). T houg h 
t he p ro blem i s s imilar to the one we have generically solved in Section [II D we cannot apply the results Eqs. ( |3.16|) , 
( |3.17 ), ( 3.18a ), ( 3.18b ), because of the limiting hypothesis (iii) mentioned in our introductory discussion Section f I D| . 
Indeed, the problem is that, in the (physically relevant) case of relativistic signals the functions a(t) and tp(t) are not 
smooth at the upper edge t = t^so- Let us see explicitly in what way they violate smoothness there. Let us first 
define, 



d_ f E'jv) 
dv \ T{v) 



(4.15) 



so that near the LSO we may write: 



E'{v) 
T{v) 



= ei(v - wlso) + 0[(v ~ i'lso) 2 ] 



(4.16) 
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If we were to use the test- mass approximation for the energy function E'(v) and the Newtonian (quadrupole) one for 
the flux function T(v) this would give 



Po(tm) 



(v) 



15 1 



2 W s (l-3i* ao )« 



27492 
4?7 



(4.17) 



We hav e num erically estimated the function ef 4 (77) = 47761(77), when using the P4-approximant of |ll| in the definition 
of Eq. (4.15). We find that to a good approximation 



4r/ef 4 (77) = ef 4 ^) ~ 26091.61194 exp (-4.474405683ri) . 



In terms of 



and using 



thso —t 

t = , r > for v < v L so , 



"lso E'(v) 1 

dv — - = -ei(w - wlso) 2 + C[(w - wlso) 3 ] 



(4.18) 
(4.19) 

(4.20) 



T{v) 2' 

and Eq. (4.2) for <f>(v), we find the following approximate representation (valid near the LSO) for the phase 4>(t): 

t ^lso = tiit , (4.21a) 



4^ 



J LSO' 



r 3/2 



(4.21b) 



Note also that Eq. ( 4.20 ) gives the following representation for v(t), and therefore for the amplitude a(t) = Cv 2 {t) 

s/2 1 



o(t) ~ a LSO 



- T 2 



2V2 1 1 

1 — T2 



We are interested in evaluating the edge contribution to the integral 



I = h+(f) 



dt a(t)e l ^f (i) = m / dr o(r) e^/ (r) . 



(4.22a) 
(4.22b) 

(4.23) 



Near the LSO boundary i.e. near the edge r = in the r-form of th e integ ral, the amplitude behaves as Eq. (4.22b) 
while the appropriate phase ipf(r) behaves, from Eqs. (4.21a) and (4.21b) as 



^+(t) ~ yj+ LSQ - 2nm(F LSO + f)r + 



4^2 



u LSO T 



where 



^/lso = ^/ (<lso) = 27r/t LS o + 0LSO ■ 



(4.24) 



(4.25) 



The appearance of fractional powers of t in the expansions Eqs. (4.22a) and ( 4.22b] ) show explicitly the violation 
of the C°° property of a(t) and ip(t) at the edge. We cannot use the integration-by-parts method to evaluate the 
expansion of / c dgc- However, we can still use the general method sketched in Section [ID. Without rotating ex plicitly 
the the r-contour in the complex plane the edge contribution to / is obtained by inserting the expansions Eqs. ( 4.22b| ) 
and (4.24) in Eq. (4.23) and expanding everything out, except for the main phase, V'/lso — 27rm(i 7 LSO + f) T which 
must be kept in the exponent. This yields 



4dge = mflLSoe^/LSO 

where, y = 27rm (F LSO + /) 



dr 1 



1 - 



2^2 1 



T2 



*e-i ^lso 



«LSO T2 



(4.26a) 
(4.26b) 
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Note that, instead of rotating r in the complex plane, we can (equival ently) consider that y possesses a small negative 
imaginary contribution: y 



y — iO. The integrals appearing in Eq. ( 4.26b| ) are evaluated by the general formula 

poo fl _j£( a +l 

y 



a + l 



-r(a + i). 



(4.27) 



This yields finally 



ma LSO e 
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LSO 



2 F LSO + / 



-iw/4 



ei VLSOy/V 



(4.28) 



The l eading contribution (oc (iy)^ 1 ) to the relativistic result Eq. (4.28) agrees with the leading contribution in 
Eq. ( 3.18b ). Note that the next-to- leading contributio n does not have the same dependence on / + Flso as the 
corre spond ing term in the non-relativistic result in Eq. ( [3. 18b ). In spite of the b reakd own of the formal expansion 
Eq. ( 3.18b ) the fractional correction given by the last term in the bracket of Eq. ( 4.28 ) is checked to be numerically 
small. This check was the main motivation for us to compute h c + sc to next-to-leading order in the relativistic case. 
The lesson is that the formal blow up of F near the LSO has only a small numerical effect on h^ 8 . This is again 
a confirmation that our results are robust under a refinement of our knowledge of the signal. We shall further check 
below that, as in the Newtonian case, has only a negligible effect on overlaps. 



C. Improved stationary phase approximation for relativistic signals 



Let us now consider the resonant contribution considered in the crucial domain where the stationary point 

is near the edge £lso As before also, the optimal approximants to h-(f) that we can construct are given by different 
analytical expressions according to the value of /. However, we need now to introduce a new definition of the two 
ranges of frequencies in which one must (minimally) divide the /-axis. More precisely, we introduce a frequency / up , 
near but below F maX j and we shall construct a "lower" approximation /j_<(/) in the range / < / up , and an "upper" 
one /i_>(/) in the range / > / up (which includes / = F max ). The optimal value of / up will be determined below. 

In the lower range, / < / up , we can draw on the work of Sec. [II D. Indeed, in that range there exists a saddle-point 
in the domain of integration. However, as that saddle-point can become rather near i max (because / up is near F max ), 
we can significantly improve the usual SPA estimate by using our previous result, i.e. by defining 



/</u P : #T(/)=C(C<(/)) 



Q ( f /) p i[lPf(tf)-Tv/4] 



Htf) 



C<(/) = -y/MW-M^) 



(4.29a) 



(4.29b) 



The label 'irspa' in Eq. (4.29a) stands for improved relativistic SPA. 

Let u s fina lly explore the optimal analytic approximation to h-(f) in the upper range / > / up . Proceeding as in 
Section IV B in this case one has 



ip f (t) ~ V> /L so + 2irm(F LSO - f)r - 



4V2 



7 lso ' 



r 3/2 



(4.30) 



where 



/LSO 



= Ipf (*LSo) = 27T/< LSO - 0LSO • 



(4.31) 



We shall use this expansion (which replaces the parabolic approximation Eq. (3.5a) used in the Newtonian case) to 
evaluate the Fourier integral Eq. ( |3.28 ). To this end we must introduce a new special function (characteristic of the 
relativistic phasing near the LSO) to replace the error function. Let us define the function 



2 ) 



(4.32) 



where the new variable f is related to r by r — ar where 
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-4/3 1/3 
LSO e l ' 



(4.33) 



In the test mass case corresponding to Eq. (4.17) the val ue of a is 49.83/(4t ? ) 1 / 3 . The value defined by the P 4 
approximant on the other hand is given by combining Eqs. ( 2.15 ) and ( 4.1 8| ) . In particular, a equals 30.055 (54.578) 
for rj — 0.25 (0.1) respectively. The index | in g$ (x ) alludes to the power r 3 / 2 replacing the power f 2 in the usual 
error function, and where the conventional coefficients 3 and 2 have been chosen to simplify some formulas (although 
they complicate others!). The final result is conveniently written in terms of a variable x given by 



2ir 

—am(F LSO - /) . 



This improved relativistic SPA is thus written as 



/ > /up : ^ S > pa (/) = ™te**t a(t LSO )gz(x) . 



(4.34) 



(4.35) 



Note that / < FLso corresponds to x > (saddle-point domain), while / > Flso corresponds to x < (absence of a 
saddle-point). Roughly speaking the variable x(f) corresponds to —((f) of the non-relativistic case, and ga (x) is the 

relativistic analogue of the combination C(()e^ appearing in the previous treatment (see, e.g., Eq. (3.31a)). 
It is useful to summarise some properties of the function g3.(x): 



5|(0) = 



1 (l-tV3). 
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x > 0, x > 1 . 
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x < , -x > 1 . 



(4.36a) 

(4.36b) 
(4.36c) 



By expanding the integrand of <?| (x) in powers of x, and integrating term by term [using the properties of the Eulcr 

T-integral after having changed the variable of integration: f = e~~z~ (u/2) § ], one proves that gi{x) is given by the 
following, everywhere convergent, Taylor-Maclaurin expansion: 



gdx) 
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(4.37) 



With about 300 terms the above series represents g^(x) accurately enough for values of x in the range x £ [—2.3, 2.3]. 
We used this series to generate the plot of <?| (x) represented in Fig. [|. Though we do not use it in this paper, note 
that for —x — > oo, the following (divergent) asymptotic expansion is also valid: 



03 (x) ~ — — > 

s w 3x ^ 

n— 



r(f ri- 



ll 



-2 



(-3x) 3 / 2 



-if(n+2) 



(4.38) 



In all our calculations of overlaps we shall define the frequency / up separating the lower range from the upper range 
by choosing a; up = 0.36 as the right hand side of Eq. (4.34). This value is chosen so that at x up one has a smooth 



transition from the lower to the upper approximation. Wc have also checked that the overlaps do not change very 
significantly for x up between 0.2 and 0.4. 

In summary our best analytic representation of time-windowed relativistic signals in the Fourier-domain would 
be defined by combining the P-approximant construction of the functions E'(v), T{v) |l5| l with the total improved 
relativistic approximants (irtot) defined as 



' irtot 



(/) - ^ kSpa + h C * 6 ° . 



(4.39) 



where hf sc is defined in Eq. ( fD| ) and /i I I spa is defined for / < / up by Eqs. (fl29a|), ( |429b| ), and, for / > / up by Eq. 
(4.35). Actually, as in the case of Newtonian signals, we have found that the inclusion of h e + 6 ° has only a minimal 



(though favorable) effect on overlaps. Moreover, such a contribution is absent in the case of real signals. Therefore, 
our final practical and best proposal consists in using only /j^ spa (For simplicity we henceforth drop the subscript 
minus). We shall henceforth refer to the improved frequency domain stationary phase P-approximants based on the 
irSPA as the SPP approximants. 
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D. Comparison between the usual SPA, the improved relativistic SPA and the 'exact' SPA (numerical DFT) 



In this Section, we test the accuracy of our analytical approximations in various ways. Fig. |^ compares an inspiral 
wave from a (20,20) Mq binary generated by three different methods: (i) directly in the time-domain and terminated 
when the instantaneous gravitational wave frequency reaches the value at the LSO (solid line), (ii) in the Fourier 
domain using the usual SPA but with a square window between / m i n =40 Hz, / max = -FLso (uSPAw) and then inverse 
Fourier transformed (dashed line) and (iii) again in the Fourier domain but using irSPA with Fourier components 
computed up to Nyquist frequency and then inverse Fourier transformed to obtain its time-domain representation 
(dotted line). We only exhibit the comparison near the crucial LSO region [Much before the LSO the uSPA is nearly 
equivalent to the irSPA and they both do a good job in representing the actual signal]. We observe that the uSPAw 
begins to get out of phase with the wave directly generated in the time-domain during the last cycle and rings a few 
times beyond the shut-off point. Our new proposal, irSPA, keeps in phase with the time-domain signal until the last 
moment although it too has a couple of low amplitude cycles beyond the LSO. 

Matched filtering involves not just the correlation of two signals but rather their weighted correlation - the weight 
coming from the detector spectral noise density. To further compare and contrast our new /-domain approximants to 
the usual frequency- windowed SPA it is conceptually useful to compare various approximations in the 'whitened-time- 



domain' introduced in Section II A above. As discussed above, in this picture (and only in this picture) the optimal 
filter consists of correlating the output of the detec tor with an exact copy of the expected signal. The whitened [i.e., 
convolved with the whitening kernel wi Eq. ( [2.10 )] signals are plotted and compared in Fig. |l^ which is the same as 



Fig. ^except that all the waves here are whitened (i.e. divided by y/ S n (f) and then inverse Fourier transformed). The 
inset in Fig. [H^ shows the full whitened signal that was originally generated in the time domain. Several observations 
are in order. First, we see how low frequency components are suppressed relative to high frequency components 
which occur in a more sensitive band of the detector. Second, we can very clearly see the non-local behaviour of 
the whitening kernel. It has the effect of softening the window imposed on the wave that was directly generated in 
the time-domain and curbing the oscillations in the irSPA beyond -Flso- Finally, this same whitening is seen to have 
worsened the mismatch of uSPAw with the whitened version of the original time-truncated signal. The conclusions 
drawn from these visual comparisons are borne out by detailed numerical experiments we performed. 



To compare the approximants more quantitatively, in Table III we list the overlaps of the exact Fourier representation 
of a model waveform (i.e., a signal generated in the time-domain and then Fourier transformed using a DFT algorithm) 
with their approximate Fourier representat ions analytically computed using one of the following: the frequency- 
windowed usual SPA [i.e., uSPAw, cf.Eq. ( |3.7| )], the im prove d Newtonian SPA [inSPA is the same as irSPAw, cf. 
Eq. ( 4.29b| )] and the improved relativistic SPA [cf. Eq. ( 4.39| )]. The uSPA and the inSPA used in computing these 



overlaps are terminated at / = F^ so , where -Fj^so is the last stable orbit frequency determined by the condition 
E' A (v) = (hence the labels SPAw and inSPAw where 'w' stands for 'windowed' — in the frequency-domain). This 

is because both uSPA and inSPA vanish at the LSO (due to the factor 1/ y - Flso) and are either not defined (in the 



case of the usual SPA) or formally vanishing (according to the definition Eq. ( 3.31a ) in the case of inSPA) beyond the 
LSO. Contrast this with the Newtonian case where it is possible to analytically extend the usual SPA beyond -Flso- 
It is generally true, as stated in DIS, that the stationary phase approximation to the Fourier transform worsens 
very significantly as we consider more massive binaries. In this sense the uSPA poorly represents the exact chirp. We 
conclude that, for massive systems with total mass m = ni\ + m,2 ^ 40M© the only uniformly acceptable analytic 
representation of the Fourier transform is the irSPA. 

V. FAITHFULNESS AND EFFECTUALNESS OF SPP APPROXIMANTS 

So far we have concentrated on developing an accurate Fourier representation of the inspiral waveform at various 
levels of approximation from Newtonian to P-approximants. In order to quantify the accuracy, we used the overlap of 
the DFT of the waveform computed using a FFT of the time-domain signal with an analytical approximation of the 



Fourier transform of the same time-domain signal using the improved SPA suggested in Sec. 111 and [V. However, an 
important question still remains: What is the total loss of accuracy due to combining the loss of precision entailed by 
the use of an analytical approximant to the FT (loss that we have shown how to minimize by defining the irspa) with 
the loss of accifracj|^] entailed by the use of some finite-order in the post-Newtonian approximation of the exact signal. 



We distinguish precision and accuracy in the same way that they are distinguished in Metrology. 
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In other words, how accurate is the approximate frequency-domain representation of a post-Newtonian approximant 
in modelling the exact FT of the exact general relativistic signal? More precisely, what fraction of the SNR of a true 
signal is the Fourier-domain approximant likely to extract? Additionally, one is also interested in knowing the biases 
induced in the estimation of parameters when using the frequency-domain approximants introduced in this work. 

We shall follow DIS in saying that a representation of a signal is faithful if it has a good overlap^ with the exact 
signal for the same values of the (dynamical) parameters (or more precisely, if the overlap is maximized for template 
parameters which have acceptably small biases with respect to the exact signal parameters). As in DIS, we employ as 
necessary criterion for faithfulness the requirement that the 'diagonal' ambiguity function be larger than 0.965. On 
the other hand, we shall say that a representation of a signal is effectual if the overlap, maximized over the template 
parameters is very near one. To use these definitions we follow DIS in introducing a fiducial exact general relativistic 
signal. In Table [V and Table |Vl| we use as fiducial exact signal the formal "test-mass case" for which the function 



E(v) is known analytically and T{v) numerically (^6). In Table [v| we use as fiducial exact signal the one defined in 
DIS for comparable masses [see Eq. (4.11) there for the definition of the exact new energy function, and Eqs. (7.1) 
(7.2) for the exact factored flux function; we took the value kq — 47/39 for the parameter defining formal higher 
PN-effects in Eq. (4.11) ]. As above we consider that the exact time-domain signal is shut off after the LSO. [For 
each considered waveform, defined by some approximate energy and flux functions Ea{v) and Ta{v), we shut it off 
at the LSO defined by the corresponding energy function Ea(v).] 

In Tables |[V] and [v| we list the overlaps for different approximants for the three 'massive' archetypal binaries 
[(1.4M Q , IOMq), (10M q , 10M Q ), and (2OM , 2OM )] that could be searched for in GEO/LIGO/VIRGO data. These 
overlaps are computed using the expected LIGO noise Eq. (1.4a) by maximising over the lag parameter r Q and 



phase (f> c but without re-adjusting the intrinsic parameters i.e., the masses of the two stars in the approximants, to 
maximise the overlap. [This implies that the value of -FLso used in the approximant is different from that in the 'exact' 
signal.] The overlaps are therefore a reflection of how accurate the various representations are in an absolute sense. In 
other words, they compare the faithfulness of the different approximants. Two independent aspects of approximation 
are investigated in these Tables. Firstly, the comparison between the two alternatives in the frequency domain: the 
usual SPA (uspaw) and our improved relativistic SPA (irspa). And secondly, the post-Newtonian order to which the 
phasing is computed. To investigate further the performance of these approximants we summarise in Table [v] the 
overlaps obtained by maximising over all the parameters in the approximants including the intrinsic ones. Thus in 
addition to maximising over the lag parameter r and the phase <f> c one also extremises over the masses of the two 
stars mi and m-i- In other words, we compare the effectualness of the various approximants. We also compute the 
bias introduced in the total mass m. 

From Tables |y|- [v| one can conclude the following: (i) The improved relativistic SPA is significantly more faithful 
and more effectual for massive systems with total mass m > 20Mq, and mandatory for m > 26Mq, (ii) Comparing 
with DIS, we see that the frequency-domain irspa does as well as the time-domain waveform even for massive binaries 
up to 40Mq; (iii) The 2.5PN SPP approximant is both a faithful and an effectual approximant for a wide range of 
binary systems (m < 40Af©). It only introduces a small bias. Note also that, in regard to effectualness, the gain in 
going from 2PN to 2.5 PN accuracy is quite significant (mainly in decreasing the biases) and especially for low-mass 
systems (which have many useful cycles), while the gain in going from 2.5PN to 3PN seems very slight. 

To summarise: // one would like to lose no more than a tenth of the events that would be observable had one known 
the exact general- relativistic signal, then the 2.5PN SPP- approximants are a must. Furthermore, unbiased parameter 
estimation requires 2.5PN SPP -approximants in all cases. 

VI. WHY ARE TIME DOMAIN RELATIVISTIC SIGNALS MORE EXPENSIVE TO COMPUTE? 

The main purpose of this work is to provide a set of tools to the experimenters so that they can generate templates 
with a minimal computational cost. We next, therefore, address the issue of computational costs of various algorithms 
for template generation. 

First, though the signal is initially given in the t ime -domain, the time-domain version of the Wiener filter contains 



a double time integration [see second form of Eq. (2J2)] which is (given the existence of FFT algorithms) much mo re 



computationally expensive than the single frequency-domain version of the Wiener filter [see first form of Eq. (2.2)] 



Therefore, in the computation of the correlation of a template with the detector output what is required is the Fourier 



11 When discussing faithfulness and effectualness we always assume, as in DIS, Eq 2.17 there, that the overlap Eq. (2.8) is 
first maximized with respect to the relative time lag (and relative phase). 
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transform of the matched filter. However, the DIS proposal was to compute the templates in the time-domain and 
compute their exact DFT using FFT algorithms. Admissibly, this procedure is still highly computation-intensive. Let 
us reason out why this is so. 

To compute the time-domain signal we need a phasing formula <f> = 4>{t). Since there is no explicit expression 
f or t h e ph asing of inspiral waves as a function of time the standard approach is to use the implicit formula, Eqs. 
(4.1 )-(4.2 ) . The binding energy E(v) and the gravitational wave flux T{v) hav e been computed, e.g. using Pade 
techniques, as explicit functions of v and these when used in Eq. (4.1) and (4.2) yield an implicit relation between 
(j) and t. However, the problem is that we need <j) at equal intervals of time (to enable us to use the standard FFT 
algorithms) and this makes the computation of <p(t) expensive: every time-sample = <f)(U) is computed by first 
solving Eq. (4.1) iteratively for Vj, the lower limit in the integral for a given ti, and then using this w, as the lower 
limit in the integral of Eq. (4.2). Though the second step is the computation of a single integral, the first step is a 
rather slowly converging (~ 10 iterations for every ti) computation. 

This problem could have been circumvented if it had been adequate to use the explicit analytical expression 
4>{t) = ba(thso — t) 5 ^ 8 + 2fc>i&fe(*LSO — t)( 4 ~ fc )/ 8 (modulo logarithms) obtained by : (i) expanding the quantity 
E'(v)/!F(v) in the integrands, in a straightforward expansion in powers of v, (ii) integrating term by term, and 
(iii) inverting analytically by successive iterations (see e.g. p4[| ). However, this straightforward PN expansion of the 
phasing formula defeats the very purpose of P-approximants and loses all the benefits brought by the constructions 
given in [Oj. Consequently, DIS had to use the iterative procedure to compute the signal phasing. By contrast, using 
(any form) of SPA, i.e. an explicit analytical f-domain expression, brings a tremendous reduction in computational 
costs. On the one hand, as we shall discuss below there is no iterative procedure involved in computing SPA. Secondly 
they are computed directly in the frequency-domain and hence lead to a further cost reduction, since time-domain 
waveforms need to be Fourier transformed using FFTs — costing N log 2 N floating point operations — in addition 
to floating point operations required to compute tim e-do main templates. 

Let us recall that the usua l SPA is given by Eq. (3.7). In this expression tf is the stationary point of the phase 
in the integral of Eq. ( 3.4b| ). At a Fourier frequency / = vj/irm the stationary point tf is given by Eq. ([L"e|), 
which is a non- iterative computation. One then substitutes this value of tf in Eq. ([0]) to compute ipf — the phase 
of the Fourier component. Moreover, the derivative of the frequency which occurs in the amplitude of the Fourier 
transform can be computed using Eq. (4.14) while the factor a(tf) oc J 2 / 3 from Eq. ( |3.2a ). Every quantity that 
appears in the SPA is computed using a straighforward integral or a mere algebraic expression. Hence, from the 
computational-cost point-of-view, it is desirable to use some SPA to generate templates. Since the usual SPA has 
been shown to be inadequate for representing time-windowed signals from massive binaries, we have proposed the 
use of corrections to (for / < / up < FLso) and analytic extensions of (for / > F up ) the usual SPA. In Table VII 
we compare for archetypal binaries, the computational costs of templates that are generated in the time-domain and 
Fourier transformed using an FFT algorithm with the computational costs for the uSPA, inSPA and irSPA. This 
Table clearly shows that it is sensible to generate templates in the Fourier domain. The SPA is up to a factor 100 
times faster and the irSPA is up t o a factor 10 times faster than the corresponding time-domain construction and 
Fourier transformation. Table VII together with Table |y| (of overlaps) demonstrates that SPP approximants while 
more expensive to generate than the usual SPA are nevertheless 'affordable', and are anyway necessary for efficient 
searches of inspiral signals in gravitational wave interferometer data. 



VII. CONCLUDING REMARKS 



After nearly two decades of detector-technology development long-baseline interferometric gravitational wave an- 
tennas LIGO/VIRGO are scheduled to become operational in about 2-4 years with target sensitivities that are good 
enough to detect inspiral events from massive (m > 20M Q ) binaries at an optimistic rate of a few per year. Searches 
are planned to be carried out over a range of 0.2-50M Q by the method of matched filtering. 

An important issue in matched filtering is the number of cycles accumulated in the correlation integral since the 
SNR grows as the square-root of the number of cycles. While this is strictly true, if the noise power spectrum of the 
instrument is independent of frequency, in practice one can only improve the SNR in proportion to the square-root 
of a "useful" number of cycles N use fa\ which is determined by a combination of the detector noise power spectrum 
and the signal's power-spectrum. We have pointed out how the number of useful cycles can be a lot smaller than the 
actual number of cycles for massive and relativistic systems: e.g. a (10M©, IOMq) [ (20Mq,20Mq)} binary system 
has only 7.6 [3.4] useful cycles in the detector's bandwidth (see Table nh. A priori , it may seem that the fewer number 
of cycles should make it easier to model the massive black hole binaries compared to the lighter neutron star-neutron 
star ones with its corresponding large number of cycles to phase. Tables |f"v|-|Vl| show that there is some truth in this, 
but that for very massive black hole binaries, these fewer cycles are in fact more difficult to model than the neutron 
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star-neutron star, or neutron star-black hole cases for two reasons: (i) they are near the end of the inspiral, i.e. when 
the radiation reaction effects drives a faster drift of the frequency which has to be modelled accurately (this is why 
we need P-approximants introduced in DIS); (ii) they might terminate due to the transition from inspiral to plunge 
while in the detector's bandwidth, and this poses the problem of accurately describing the Fourier transform of a 
time- windowed signal (this requires the correction factors introduced in this paper). All this places stringent demands 
in modelling the waveform in the Fourier-domain and due attention needs to be paid to delicate issues of detail. This 
task is all the more important that the first detections expected from LIGO/VIRGO are likely to concern massive 
systems with m ~ 25 ± 5Af Q , for which the LSO frequency lies near the middle of the sensitivity curve [see Fig. 
To this end, the present work makes two new robust (i.e. assumption-independent) contributions: 

• the proposal of stationary phase P-approximants (SPP) which combine the excellent performance of our time- 
domain P-approximants Jl3[ with the analytic convenience of the stationary phase approximation without serious 
loss of event-rate. These Fourier-domain P-approximants perform as well as their time-domain counterparts in 
extracting the true general relativistic signal. 

• the definition of a universal Newtonian-like 'edge-correction' factor C(£(/)), as well as its relativistic complement 
9z(z(f)) which take into account the frequency-domain effects, concentrated around (and on both sides) of 
^max = f (imax)) for signals which are abruptly shut off, in the time-domain, after t max . 

In addition to these new achievements, let us mention two other useful contributions, of a more technical nature: 
(i) our recommendation to systematically use a smooth time window at the lower frequency side to conveniently 
and efficiently suppress spurious oscillations due to a numerical low frequency cutoff and (ii) the emphasis on the 
comparison of the form of the signals in the 'whitened' time-domain. 

Based on the detailed analysis presented in this paper we find that for post-Newtonian template generation of binary 
systems of total mass m < 5M it suffices to use the usual SPA (without correction factor) of the P-approximants 
defined in DIS. On the other hand, in the total mass range 5M Q < m < 40M Q , it is crucial to use our new SPP 
approximants to construct the frequency-domain templates. 

In addition to the construction of the SPP approximants, the paper has examined in detail the Fourier-domain 
effects entailed by a sharp time-domain windowing. As emphasized in the introduction, at our present stage of 
knowledge, one cannot be sure that a template waveform terminated (in the time-domain) at the LSO is an accurate- 
enough representation of a real GW signal coming from massive binaries (say with m < 40M Q ). We have given several 
plausibility arguments towards justifying this assumption: brevity of the plunge, and an expected frequency separation 
from the merger signal. In the absence of knowledge of the transient plunge signal and of the final merger signal, we 
have argued that it is best to use a template waveform which is terminated at the LSO. [Actually, we anticipate that 
the effcctualness of the template waveform will be increased if we allow it to be terminated at a frequency somewhat 
larger than -Flso (thereby allowing it to approximately represent the plunge waveform) .] Consequently, this work has 
concentrated on signal models that are truncated in the time-domain by a step- function and has aimed at constructing 
the best associated Fourier-domain analytical representation for this possibility. 

We have also pointed out that the opposite assumption of an abrupt termination at Flso of the usual SPA in the 
frequency domain implies, when viewed from either the time-domain or the whitened time-domain, the existence of 
some coherent oscillations 'ringing' after the LSO crossing. We have done another numerical experiment on this issue, 
by appending to the inspiral signal a smooth decay taking place over less than <iF^g time-scales. We have found 
that our improved SPA was a reasonably good representation of (the FT of) such a signal, and definitely a better 
representation than the usual SPA one. Let us finally reiterate that, we do not claim to have conclusively ruled out 
the possibility that a frequency-windowed SPA may perform better compared to the time-windowed SPA we propose 
here. This important issue is not settled though we conjecture that this is unlikely. Anyway, this paper is the first 
one to explicitly construct the frequency domain version of the time-domain P-approximants which were shown in 
DIS to bring indispensible improvements over the usually considered T-approximants. Therefore, even in the unlikely 
case where a straightforward frequency-window turns out to be a better model than the time-window assumed in 
most of this work, one will still require the formulas given in this paper (with the trivial change of replacing the 
correction factors C(Q by a 9 function 0(Flso — /)) to generate sufficiently accurate f-domain filters. In view of these 
comments, we feel there is a urgent need to model more precisely the transition from the inspiral signal to the plunge 
signal 0] close to the last stable orbit. We hope that the techniques (if not all the details of our construction) used 
in this work to handle the blow up of F(t) at the LSO will be useful (maybe with some modifications) even if, on a 
later examination, this blow up turns out to be an artefact of an approximation which may drastically alter with a 
better treatment of the transition to the plunge. Only with this improved understanding and its implications for the 
construction of templates can one build even more optimal templates for massive binaries and maximise our chance 
of detecting them. Independently of issues such as windowing in time versus windowing in frequency or the nature of 
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the plunge we feel that in general P-approximants are much better tools than the Taylor approximants. We hope to 
come back to this question in a future work . 

Another aspect that needs to be looked into is the issue of whether whether the interferometers will work in the 
time-domain or the frequency-domain. If indeed, they would decide to work in the time domain: i.e., to store 
the raw output, and to transform it nearly online in the dcfiltcrcd time-domain equivalent GW amplitude h(t) the 
analysis of this paper would be irrelevant. In that case, one should store the Wiener transformed time-domain filter 
K(t) = wi * h(t). However, with the presently available computational resources it seems hopeless to filter in the 
time-domain. We therefore anticipate that, though the raw detector output will be stored in the time-domain, all 
filtering will be done in the Fourier domain. In this event, the robust aspects of the present analysis will be relevant 
even if not the details. 

The formalism developed in this paper can be applied not only to initial interferometers but also to future generations 
of interferometers. We have refrained from applying our formalism to the case of LIGO II since the LIGO II design 
is at the moment in a state of flux and any quantitative results we may quote will soon be irrelevant. However, we 
should expect the results of this work to be important for any detector that works with a lower seismic cutoff and a 
broader bandwidth than LIGO I, since in such cases we will have to match the signal's phase for a larger number of 
effective cycles. 

There are several notable and obvious improvements that need to be pursued. The sensitivity to the value of 
-Flso needs to be investigated [in particular, our improved SPA will probably maximize their overlaps with the real 
signals if we allow some flexibility in the choice of FLso (within some limits)]. Once the results of 3PN generation 
of gravitational waves are available Q and are combined with the 3PN results on the dynamics Q they must be 
included in the construction of templates. In our discussion we have not considered waveforms from binaries with 
spinning compact objects. Nor have we included the effect of eccentricity |30| on the detectability p3[ . These are 
unarguably important physical effects that need to be incorporated in later data analysis algorithms. Future research 
in this area should shed light on these issues. 
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APPENDIX A: LIST OF SYMBOLS 



a(t) 
a 

cspa 

C n {t\ — t 2 ) 

c(0 
s 

V 

erfc(x) 

E(v) 
£l 

£2 

FFT 

f-domain 

f-window 

F(t) 

Hv) 

^rnin (-^max ) 

^LSO 

-^Nyquist 



GW amplitude; h(t) = 2a(t) cos <f>{t) 



-4/3 1/3 
2 U LSO e l 



Eq. O 



corrected SPA; Eqs. (|3T29a| )- ([3.29b|) 
correlation function of noise 

ierfc(e l7r / 4 £); correction factor; softened step function 



leading phase correction to SPA; Eq. (3.9) 
symmetric mass ratio = mii7i2/(mi + 1TI 2) 2 
complementary error function; Eq. ( 3.21 ) 



A. ( E '( v 

dv y 3~{v, 



Eq. O 



dimensionless energy function 
= Eq. (P) 



1 a(t)<t>(t 

= liWj - I 1 Ht) I = 1 . 
— I 02(^)1 I 2ir F 2 (t) I 2ttN' 

Fast Fourier transform 
frequency-domain 
frequency window 
instantaneous GW frequency 
flux function 

GW frequency at t min (t m!iX ) 
GW frequency at LSO 
Nyquist Frequency 



Eq. (3.8) 
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/det 
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/up 

9% ( x ) 
GW 

r 

W) 

K(f) 
h(t) 

kf) 

M/) 
h-U) 

/i inspa (/) 

^ + dSC (/) 
h intot (f) 
h}l Spa {f) 
h Mot (f) 

inspa 

irspa 

l(v) 

rn 

M 

n(t) 

JVtat 

AT(F) 

^now(F) 
■^"useful 

o 
p 

SPA 

Sn{!) 

a(t,zi,z 2 ) 

^min 
^max 
T n 
T 

6 

uspa 

uspaw 

uspan 

v 

VM 
*»(/) 

toi(t) 
ttii (r) 

2 v ' 

x 



Fourier frequency 

characteristic detection frequency; minima of effective GW noise \J fS„ (f) 
frequency at which e?SNR /<i(ln /) peaks 
seismic frequency 

transition frequency between the low and high frequency approximations for the irSPA 
= f™dre< 3xf - 2f3/2 ); Eq. (p2|) 
Gravitational wave 
Gamma function 

squared amplitude of effective GW signal; = N(f)a 2 (f) 
squared amplitude of effective GW noise; = fS n (f) 
time domain signal 

Fourier transform of h(t); h(f) = dte 27Tift h(t) 
Fourier transform of non-resonant part of h(t) 
Fourier transform of resonant part of h(t) 

improved Newtonian SPA corresponding to h-(f); Eqs. ( 3.27 ) ( [3. 31a ) 



edge approximation to h+(f); Eq. ( |3.18b| ) 

= h e _^ se + /i™ spa : total improved Newtonian SPA of h(t) 

improved relativistic SPA corresponding to h-(f); Eqs. ( |f.29a|) , ( 4.29b ), (|4.35|) 

= h+ dgc + /i i ! spa : total relativistic SPA of h(t) 

improved Newtonian SPA 

improved relativistic SPA 

factored flux function 

total mass of the binary 

chirp mass = rf/^m 

noise 

whitened noise; = w i (t) * n(t) 
total number of cycles; Eq. ( 2.18 ) 
instantaneous number of cycles; Eq. ( 2.1 9| ) 
instantaneous number of cycles in Newtonian case; Eq. (|2.26| ) 
instantaneous number of cycl es in relativistic case; Eq. (2.28) 
useful number of cycles; Eq. (|2.24|) 



overlap (Normalised Ambiguity Function); Eq. (2.8) 
GW phase 

P-approximant of order v n 
signal to noise ratio 
stationary phase approximation 
two-sided noise power spect ral d ensity 
weight function in p 2 ; Eq. ( 2.34 ) 
smoothing time window; Eq. (3.35) 
starting time of the signal 

time at which the signal terminates or is terminated 
Taylor approximant of order v n 

tLSQ—t 
m 

Heaviside step function 
usual SPA 

usual SPA frequency windowed 
usual SPA up to Nyquist frequency 
invariant velocity (mnF) 1 ^ 3 
invariant velocity (ttMF) 1 / 3 
weight factor °^ W: 

correlation inverse of n oise correlation function; Eq. ( |2.3| ) 
whitening kernel; Eq. (2. id) 



Coif) 



f am( F LSO -f); Eq. (^34| 



irF{t f )(t f -t max ); Eq. (3.29b) 
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C<(/) = -y/MW-Mtm^); Eq. §3 

c>(/) = ^= 1 ; Eq. (HI) 
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TABLE I. Number of useful cycles for some representative systems and approximants. For each of the cases the corresponding 
total number of cycles is listed within rounded brackets in the line below. 



System Newtonian Relativistic:P4 

1.4 Mq - 1.4 M Q 173 169 

(1588) (1586) 
1.4M -1OM Q 51 37 

(348) (320) 
10 Mq - 10 Mq 12 7.6 

(57) (47) 
20 Mq - 20 Mq 9.2 3.4 

(15) (10) 



TABLE II. Accuracy of the various approximations to the Fourier transform of the chirp signal in the test mass limit for 



Newtonian waveforms. The accuracy is estimated via the overlap, Eq. (2.8), of a waveform generated in the time-domain 
and then Fourier transformed using an FFT algorithm, with waveforms of exactly the same parameters generated directly 
in the frequency domain via one of the following approximations: The usual stationary phase approximation truncated at 
/max = fkyquist, (uSPAn, column 3), the usual stationary phase approximation truncated at /max = Jlso (uSPAw, column 
4), the corrected SPA (cSPA, column 5) and the improved SPA (inSPA, column 6), both extended up to /Nyquist, for different 
values of the total mass (column 1). The corresponding Jlso is listed in column 2. [We deal only with equal mass binaries 
r\ — 1/4.] In the last column we also list the total overlap 'intot' to exhibit the relative importance of the 'non-resonant' 
contribution. In this table and all others the smooth time-window starts at 30 Hz; the seismic cutoff for LIGO is at / s = 40 
Hz. 



tji/Mq 


-Flso/Hz 


</£ FT ,^ SPAn ) 


(^ FT ,feN SPAW ) 




(^ FT ,^ nSPA ) 




70.0 


63 


0.1536 


0.6361 


0.9231 


0.9763 


0.9944 


60.0 


73 


0.2294 


0.7302 


0.9489 


0.9721 


0.9873 


50.0 


88 


0.3336 


0.8062 


0.9724 


0.9824 


0.9934 


40.0 


110 


0.4708 


0.8589 


0.9862 


0.9952 


0.9993 


30.0 


147 


0.6682 


0.9214 


0.9964 


0.9987 


0.9977 


20.0 


220 


0.8811 


0.9681 


0.9968 


0.9974 


0.9986 


15.0 


293 


0.9431 


0.9838 


0.9999 


0.9998 


0.9997 


14.0 


314 


0.9528 


0.9868 


0.9995 


0.9997 


0.9993 


13.0 


338 


0.9641 


0.9900 


0.9986 


0.9989 


0.9993 


12.0 


366 


0.9700 


0.9912 


0.9996 


0.9997 


0.9995 


10.0 


440 


0.9839 


0.9953 


0.9989 


0.9990 


0.9994 


5.0 


880 


0.9983 


0.9988 


0.9991 


0.9991 


0.9992 


3.0 


1466 


0.9987 


0.9987 


0.9985 


0.9985 


0.9984 
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TABLE III. Accuracy of the various approximations to the Fourier transform of a chirp signal computed in the test mass 
limit for P-approximants at different post-Newtonian orders v n (column 1) is measured by their overlaps with a waveform 
of exactly the same parameters but computed in the time-domain and then Fourier transformed using an FFT algorithm. 
The approximations considered are the usual SPAw (column 2), inSPAw (column 3) and irSPA (column 4) (with the choice 
x up ~ 0.36). Both uSPA and inSPA are windowed in frequency beyond Flso- Only the irSPA has a spectrum extending beyond 
^lso though we put a numerical cutoff at Kcutoff = —20. x up — 0.36 and x cuto g = —20 in all subsequent tables. 



n 


uSPAw 


inSPAw 


irSPA 










mi = 1.4M©, m 2 = 1OM 


4 


0.9967 


0.9986 


0.9994 


5 


0.9965 


0.9990 


0.9997 


6 


0.9965 


0.9986 


0.9993 


mi = m,2 — IOMq 


4 


0.9764 


0.9762 


0.9951 


5 


0.9771 


0.9775 


0.9955 


6 


0.9751 


0.9786 


0.9953 


mi = m 2 = 20M© 


4 


0.8613 


0.8613 


0.9891 


5 


0.8680 


0.8773 


0.9819 


6 


0.8897 


0.9038 


0.9829 



TABLE IV. Faithfulness of the Fourier domain P-approximants in the formal test mass case. Values quoted are the minimax 
overlap (see DIS) of an approximant waveform generated directly in the frequency-domain with the exact waveform generated 
in the time-domain and Fourier transformed using the FFT algorithm. The approximations considered are the usual uSPAw 
(column 2) and irSPA (column 3) As mentioned earlier x up = 0.36, x cuto B = —20. For compactness of presentation, in this 
Table and Tables [v| and [vi] we use for the column headings the abbreviation uSPAw to denote (hy? T , /tp^, PAw ) irSPA to 
denote (/£ FT ,/^f A ). 



n (1-4,10) 


(10,10) 


(13,13) 


(20,20) 


uSPAw irSPA 


uSPAw irSPA 


uSPAw irSPA 


uSPAw irSPA 


4 0.8360 0.8316 

5 0.9755 0.9729 

6 0.9921 0.9903 


0.9596 0.9707 
0.9727 0.9914 
0.9739 0.9938 


0.9513 0.9965 
0.9473 0.9973 
0.9479 0.9972 


0.8248 0.9785 
0.8283 0.9855 
0.8285 0.9856 


TABLE V. Faithfulness, as in Table IV but for waveforms constructed from post-Newtonian formulas for a binary with stars 
of comparable masses. The value of the parameter ko = 47/39. 


n (1-4,10) 


(10,10) 


(13,13) 


(20,20) 


uSPAw irSPA 


uSPAw irSPA 


uSPAw irSPA 


uSPAw irSPA 


4 0.7919 0.7898 

5 0.9765 0.9736 

6 0.9965 0.9958 


0.9596 0.9584 
0.9820 0.9924 
0.9835 0.9947 


0.9485 0.9665 
0.9657 0.9921 
0.9669 0.9921 


0.8764 0.9790 
0.8831 0.9815 
0.8860 0.9868 
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TABLE VI. Effectualness of the Fourier domain P-approximants in the formal test mass case. We list minimax overlaps of 
the approximate waveforms generated directly in the frequency-domain with the exact waveform generated in the time-domain 
and Fourier transformed using an FFT algorithm. However, in addition to maximisation over the lag parameter r and the 
initial phase of the approximate waveform we also maximize over the two masses mi and m-2. The percentage bias in the 
estimation of the total mass 100(1 — mA/m) are listed in parenthesis below the overlaps. 



n (1-4,10) 


(10,10) 


(13,13) 


(20,20) 


uSPAw irSPA 


uSPAw 


irSPA 


uSPAw irSPA 


uSPAw 


irSPA 


4 0.9942 0.9963 
(-2.300) (-2.744) 

5 0.9882 0.9987 
(+0.480) (-1.301) 

6 0.9970 0.9982 
(-0.311) (-1.166) 


0.9798 
(-0.664) 

0.9795 
(-0.626) 

0.9805 

(0.000) 


0.9984 
(-0.858) 

0.9979 
(-0.467) 

0.9964 
(-0.626) 


0.9651 0.9979 
(-0.962) (-1.081) 

0.9653 0.9968 
(-0.631) (-0.481) 

0.9649 0.9963 
(-0.662) (0.000) 


0.9166 
(-2.500) 

0.9174 
(-1.875) 

0.9174 
(-1.250) 


0.9966 
(-1.250) 

0.9965 
(-0.340) 

0.9966 
(-0.340) 


TABLE VII. In this Table we list times required to generate templates first in the time-domain and then Fourier transforming 
using an FFT (column 2) and compare them with times required to construct the same templates directly in the Fourier domain 
using one of the three approximation schemes: uSPA (column 3), inSPA (column 4) and irSPA (5). For each post-Newtonian 
family (column 1) the time tp^ T required to compute time-domain P-approximant templates is the highest and we normalise 
all times by this value. 


n FFT 




uSPA 


inSPA 




irSPA 






r n. I r n 


*p„ SPA Ap! T 




SPA nkVl' 

n 1 Pti 






mi = 


= 1.4M , m 2 = 1OM 






4 1 
6 1 




0.013 
0.014 


0.021 
0.021 




0.089 
0.090 


mi = m,2 = 10M Q 


4 1 
6 1 




0.009 
0.009 


0.015 
0.016 




0.11 
0.11 
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FIG. 1. This Figure compares the signal-to-noise rati o (S NR) for an inspira l sign al searched by means of two different 
filters, in two different interferometers [initial LIGO, Eq. (1.4e) or VIRGO, Eq. (1.5a)]. The exact signal h is assumed to be 
a time-truncated Newtonian chirp (sufficiently well approximated by the inSPA). The solid lines represent the optimal SNR, 
obtained when the filter k is identical with the exact signal h. The dotted lines represent the sub-optimal SNR, obtained when 
the filter k is the frequency- windowed usual SPA (uSPAw). In the case of LIGO, for low mass binaries uSPAw extracts the 
full SNR, but at higher masses, which are the most crucial ones for initial interferometers, it loses SNR up to a factor of 1.5 
leading to a loss in the number of detectable events up to 70%. In the case of VIRGO, uSPAw performs quite well when the 
total mass m < 50Mq, and requires inSPA for heavier binaries. The loss in the number of events caused by uSPAw for heavier 
mass binaries in the range m > 50Mq is unacceptable. 
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FIG. 2. The whitening kernel, Wi (t), Eq. (2.1C ), is shown [for the initial LIGO noise curve, Eq. ( |l.4a| )] plotted as a function 
of time. We see that it is quite sharply peaked at t = 0. The peak drops quite rapidly as we move away from the origin thus 
indicating that we have almost local-in-time filters. The curve also indicates the importance of knowing the plunge signal over 
a time-scale of several 10's of milliseconds so that the inspiral signal can have a good overlap. 
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FIG. 3. The Newtonian instantaneous number of cycles N(f), the square of the amplitude o 2 (/), their product N(f)a 2 (f), 
the reciprocal of the effective GW noise = fS n (f), and d(S N R) 2 / d(log f) are plotted as a function of /. The scale on the 

y-axis corresponds to N(f) on the left and all other functions are plotted on an arbitrary scale indicated on the right. The top 
panel is for a lighter mass binary (mi = IAMq, 7712 = 1OM ) and the lower one for a heavier mass system (mi = rri2 = 1OM ). 
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FIG. 4. Instantaneous number of cycles, Eq.( p.l9| ), near the LSO as a function of time for the Newtonian case, Eq. (2.26), 
and the relativistic cases, Eq. (2.28), defined by the 2PN approximant P4. Also plotted is the development of the wave- form 
hp i (t) on an arbitrary scale. The plots show how the number of useful cycles diminishes as one gets close to the LSO. The 
Figure refers to mi = 7712 = 2QMq. 
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FIG. 5. The real and imaginary parts of the complex correction factor C(() — in terms of which the new frequency domain 
approximants are represented in the Newtonian-like cases — as a function of £. The real part is a 'softened step' while the 
imaginary part an oscillatory function vanishing at the origin and at large positive and negative values of its argument. 




45 



FIG. 6. The power per logarithmic bin of the squared SNR d(p 2 ) /d(log /) = f\h(f)\ 2 /S n (f) for an arbitrarily normalised 
Newtonian signal computed from its DFT and its various approximate representations computed up to the Nyquist frequency: 
uSPAn, cSPAn and inSPAn. In the most sensitive range of frequencies, our final proposal inSPAn agrees with the FFT quite 
well. The uSPAn grossly overestimates the true signal power at frequencies / > Flso- The last stable orbit frequency in 
this case (mi = mi = 10M Q ) is at about 220 Hz (the vertical line). Observe that, therefore, even uSPAw would significantly 
overestimate the signal power up to Jlso- 
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FIG. 7. The instantaneous GW frequency F(t) versus t for the Newtonian and relativistic cases during the last few orbits 
before the plunge at LSO. Notice the rapid increase in the inspiral rate close to the last stable orbit in the relativistic case. 




FIG. 8. The real and imaginary parts of the <?3 (x) function Eq.(4.37) in terms of which the improved relativistic SPA 
(irSPA) is represented. The thick lines represent the asymptotic behaviour as x — > ±00 given by Eqs.( 4.36b) and (4.36c). 
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FIG. 9. Visual comparison of a chirp generated directly in the time-domain (solid line) with the inverse Fourier transforms 
of the usual SPA terminated at LSO (uSPAw) (dashed line) and the improved relativistic SPA (irSPA) extended up to the 
Nyquist frequency (dotted line), all during the last few cycles before the last stable circular orbit frequency is reached. The 
Figure corresponds to the relativistic (second post-Newtonian P-approximant case (P4)) inspiral of a binary of total mass 
m = 40M©. 
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FIG. 10. This plot is the same as in Fig. |9| except that we compare whitened signals to show the effect of the detector 
response function on the time development. The inset shows the entire time-domain signal starting from Fgw = 30 Hz and 
terminating at FLso- 
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